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FOREWORD 


This final report, prepared by the Lockheed Missiles and Space Company for 
the National Aeronautics and Space Administration-Lewis Research Center, 
documents the technical effort for Contract NAS3-17798* "Numerical Simulation 
of Low Gravity Draining." (The Contract started on Oct. 9 , 1973 and was 
completed on May 10, 1976. ) Mr. E. P. Symons was the NASA Project Manager. 

The authors wish to thank Mr. Symons for support throughout the study and, 
in particular, for supplying Lewis Research Center experimental data for 
comparison with our analytical and numerical results. 

The present study is a continuation of an earlier study, "Low Gravity 
Draining from Heraispherically Bottomed Cylindrical Tanks, " performed by 
P. Concus, H. M. Satterlee and the authors under Contract NAS3-11526. The 
problem formulation section of that study, reported in reference 1, is 
essentially reproduced herein. The authors acknowledge, therefore, the 
previous contributions of Drs. Concus and Satterlee . < 
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ABSTRACT 


Assuming that the liquid is inviscid and incompressible, that its motion is 
irrotational and axisymmetric, and that its solid- liquid contact angle is 
constant (5°), the initial boundary -value problem is solved numerically. 
Excessive mesh distortion encountered with strictly Lagrangian or Eulerian 
kinematics is avoided by introducing an auxiliary kinematic velocity field 
along the free surface to vary the trajectories used in integrating the 
ordinary differential equations simulating the moving boundary. The compu- 
tation of the velocity potential ij based on a nonuniform triangular mesh, 
which is automatically revised to varying depths to accommodate the motion 
of the free surface. These devices permit calculating draining induced 
axisymmetric slosh through the many (or fractional) finite amplitude oscil- 
lations that may occur depending upon the balance of draining, gravitational, 
and surface tension forces. Velocity fields, evolution of the free surface 
with time, and liquid residual volumes are computed for three and one half 
decades of Weber number and for two Bond numbers, tank fill levels, and 
drain radii. Comparisons with experimental data are very satisfactory. 
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SUMMARY 


The problem is to simulate low gravity draining from a hemispherically bot- 
tomed tank over as wide a range of parameters as possible and to compute the 
residual volume at vapor ingestion (the entrance of the free surface into the 
drain). Fourteen simulations are presented covering the contractual nondimen- 
sional data ranges: two drain radii, l/30 and l/lO of a tank radius; two 

initial volumes, 5 n f 3 and 8m/3; two Bond numbers (ratio of gravitational to 
surface tension forces), B=0 and 5; and three and one half decades in Weber 
number (ratio of draining to surface tension forces), W=.01 through 50. The 
range in Weber number extends that of experiment by two decades. Simulation 
agrees very well with experiment; application of qualitative results from this 
study explains both areas of excellent agreement and regions of divergence in 
two comparisons with experiments performed at NASA-Lewis Research Center. 

Results include plots, in a form suitable for use in design studies, of residual 
volume and vapor ingestion time against W and W/(l+B), a simple correlation 
parameter occurring in the transformation between nondimensional time scales for 
draining and sloshing problems. The S-shape of the residual volume plots is 
indicative of three draining regimes: (i) a capillary regime with low Weber 

numbers, low residuals, and many waves having small amplitude; (ii) a draining 
dominated regime with high Weber numbers, large residuals, and a fractional 
wave in which much liquid is dynamically trapped on the wall as the center is 
rapidly withdrawn; and (iii) a transition regime with few waves of large ampli- 
tude which may become very irregular in zero gravity. Increasing the Bond 
number may (i) increase the number of waves and decrease their amplitude, 

(ii) change the motion from draining dominated to oscillatory, or (iii) reduce 
the amount of liquid trapped on the wall. Drain radius affects residual volume 
for W < < 1 or, at a fixed W, for B > > 1. Increasing initial volume when wave 
motion occurs gives more waves of the same nature; however, increasing it in the 
draining dominated regime gives gravitational and capillary forces more time, 
before vapor ingestion, to drain the liquid trapped against the wall into the 
bulk flow. In a simulation for B=5, W=10, and initial volume 8 m/3, the liquid 
trapped on the wall evolved into a thin wall sheet with the flow rate down the 
wall less than that in the bulk flow. The program terminated because the separa- 
tion into two flows had not been anticipated. There may be problems near the 
boundary of the draining regime for which the validity of the inviscid model may 
be questioned. However, in a survey program it is appropriate. 
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The model is based upon the irrotational flow of an incompressible, inviscid 
fluid. The model assumes that the liquid is initially at rest, includes the 
static and dynamic effects of surface tension with the solid-liquid contact 
angle taken as constant (here 5°)» retains the nonlinear convective terms in 
the differential equations derived from the Bernoulli equation, and assumes 
that, after impulsive opening of the drain, the outflow rate is uniform and 
constant across the orifice. The mixed initial boundary-value problem result- 
ing from the nondimensionalization and specialization of the model to axi- 
symmetric motion in a hemispheric ally bottomed circular cylindrical tank is 
reformulated as an initial value problem for three ordinary differential equa- 
tions for the coordinates and velocity potential at a point on the moving 
boundary with the connection to the drain given by the solution of Laplace's 
equation within the liquid. The reformulation shows that almost all the 
action in the problem is at the free surface, a fact which motivates the method 
of numerical solution and the form of the graphical presentation of the data. 

A significant result is the demonstration that the method of moving triangular 
meshes supplemented by the method of variable trajectories, is a practical and 
economic way to solve thi6 difficult moving boundary problem. A finite differ- 
ence analog of Laplace's equation is solved on a nonuniform triangular mesh 
which is revised as needed, deeper and deeper into the liquid, in response to 
the changing shape of the free surface. Mesh rows are dropped to accommodate 
outflow. Parts of mesh columns are dropped or added in response to changing 
thickness of the wall sheet. The method of variable trajectories was developed 
specifically for this study to deal with the mesh difficulties encountered as 
the free surface stretches and contracts when either Lagrangian or Eulerian 
kinematics are used. A 6mall auxiliary kinematic velocity directed along the 
free surface and chosen to maintain a desired mesh distribution is added to 
the Lagrangian velocity components at each time step; the differential equation 
for the velocity potential is appropriately modified to account for the change 
in direction. 
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INTRODUCTION 


This analytical and numerical study of low gravity draining from hemi- 
spherically bottomed, cylindrical tanks is a continuation of an earlier 
study [l]. The mathematical formulation of the problem, the basic ap- 
proach to its numerical solution, and much of the computer programming are 
carried over from the previous to the present work. Successful simula- 
tions were obtained in the former for small and moderate fill levels for 
Weber number one and greater, which, at that time> corresponded to values 
obtainable in the short draining time available within the period of 
weightlessness provided by the 2.2- and the 5- to 10-second Zero-Gravity 
Facility at NASA Lewis Research Center [2-5]. Comparisons of the computed 
and experimental data available led to the conclusion that a model based 
upon the irrotational motion of an incompressible, inviscid fluid satis- 
factorily describes the fluid dynamics of low gravity draining. The model 
assumes the liquid is initially at rest; includes the static and dynamic 
effects of surface tension with the solid-liquid contact angle taken as 
constant (here 5°); and retains the nonlinear convective terms. The drain 
rate is taken as constant and uniform across the orifice. The mixed 
initial-boundary value problem resulting from the specialization to axi- 
symraetric motion in a hemispherically bottomed, cylindrical tank is de- 
tailed in the next section, which is essentially the nondimensional prob- 
lem formulation of the previous report repeated for convenience and com- 
pleteness (cf. Figure l). 

The fundamental goal of the present study was to simulate low gravity drain- 
ing over as wide a range of parameters as possible to determine residual 
volume as a function of Weber number. Bond number, fill level, and drain 
radius. Particular emphasis was directed toward developing the program of 
the previous study in three areas: low Weber number cases, where draining 

induced' wave motion predominates; high Weber number cases, where develop- 
ment of a sheet along the wall is the dominant feature; and high fill level 
cases, which require long-time simulation. Only those changes in the 
numerical scheme or its implementation that directly contribute to achiev- 
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ing this goal have been made. Details are given in the subsequent section. 
Numerical Analysis; however, some are sketched here. 

Both Eulerian and Lagrangian kinematics on the free surface failed in -the 
previous study by producing large mesh spacings near the wall when a wall 
sheet developed; and the latter also failed by pushing me6h points too near 
the axis, thereby exciting instabilities inherent in the singularity of the 
coordinate system along the axis. 'Do eliminate these difficulties, an 
auxiliary kinematic velocity directed along the free surface is used to vary 
the trajectories of the mesh points on the free surface so that the ratio 
between adjacent mesh intervals remains nearly constant. Ihe Bernoulli . 
equation is modified to accommodate the change in the direction of inte- 
gration. This method of variable trajectories includes both the Eulerian 
and the Lagrangian formulation of the problem. The latter was introduced 
for low gravity reorientation by Moore and Perko [6]. 

The method of moving triangular meshes, carried over from the previous 
program, is now capable of carrying a low Weber number problem through 
many waves and of generating satisfactory meshes for more than 78 OO time 
steps (cf. Figures 2 and 3 )* It also handles the development and thinning 
of a wall sheet adequately. The method of truncated direct solution for 
determining the velocity potential and the representation of the free sur- 
face by cubic spline interpolation are carried over from the previous work. 

The behavior of the streamlines with time reflects the fact that almost all 
the action is at the free surface. Well up into the tank, they are fixed; 
as they approach the free surface, they begin to vibrate; and they reach 
their maximum amplitude at the free surface. Thus the trajectories of the 
points in which the free surface meets a fixed set of streamlines is a 
natural way to present their behavior with time. Moreover, these points are 
easily determined in the course of computing the normal derivative. 

Computer plots are a large part of the report (cf. Figures 4-17) : for each 
case there is a figure with the time history of the free surface and the 
associated velocity field presented on facing pages. On the first page 
successive positions of the free surface are shown. The eleven lines cross- 
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ing the free surfaces are the trajectories of the points in which the free 
surface meets a fixed set of streamlines (.02, .1, .2, . .., .8, .9, .98). 

The reader is cautioned that these trajectories are not streamlines. The 
trajectories do reveal the changes in motion that are taking place under- 
neath the rather smoothly falling free surface. When compared case to case, 
they show how features of the low gravity draining phenomenon vary with 
parameter changes. Moreover, the trajectories can be used, with caution, to 
make statements about the nature of the flow. For example, the central tra- 
jectory divides the free surface into two parts, a central disk and an outer 
annulus, each of which furnishes one half of the flow out the drain. When 
the central trajectory is far to the right of its nominal position (r=.707 . ..), 
the average velocity in the outer annulus must exceed that in the central 
disk. Because the draining free surface is generally lower at the center 
than at the wall, when a streamline end is moving to the right, a velocity 
component is probably developing down the free surface; when it is moving to 
the left, the component along the free surface is being choked off. Examples 
of streamline plots are included in many figures to illustrate specific 
points and to provide the reader with a way to check his use of these tra- 
jectories. Time histories of the centerline and wall points on the free sur- 
face are given on the second page of each figure. A free surface of interest 
can be associated with a time (or vice versa) by using the height at the 
centerline or wall as a common entry in the graphs on the facing pages. 

The simulations complement the' experimental work in low gravity draining at 
Lewis Research Center [2-5]. Agreement between experiment and simulation is 
remarkably good in one case with both parameters and initial conditions nearly 
identical (cf. Figure 18) and is acceptable in another case even though 
the experimental measurements at the wall are in doubt (cf. Figure 19). 

The simulations seem to have shorter draining times and, therefore, higher 
residuals than expected from experiment. Two explanations are proposed. The 
experimental technique may produce an initial velocity distribution that gives 
longer draining times than the zero initial velocity distribution of the simu- 
lations. The five degree contact angle condition assumed in the simulation 
may retain enough liquid high up near the wall to account for the discrepancy. 
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Both hypotheses can be tested numerically; it is possible to materially re- 
duce the contact angle and to simulate one or more full drop tower runs. 

Five of the fourteen simulations are for Weber numbers .01 and .1 which ex- 
tend data two orders of magnitude below the experimental range with conserva- 
tive residuals. 

It would be surprising if some evidence of waves did not appear on the wall 

sheets; Figure 13 (c) and (d) shows such a case. In general these are sup- 

pressed by the large mesh spacing used on the free surface. In one long- 
time simulation a wall sheet develops, begins to fall, thins out, and then 
begins to neck off when it is nearly a tank radius long and several hundredths 
of a radius thick. Figure 16 shows this development in detail and the final 
pages of the figure show a remarkable example of automatic generation of a 
mesh. This suggests that there may be combinations of Bond number, Weber 
number, and fill level for which liquid may be trapped on the wall. Further 
investigation is needed. The fourteen simulations provide information about 
the fluid dynamics that is difficult to measure or observe in tests and 
thereby give a basis for developing a qualitative explanation for low gravity 
draining phenomena. 

Computed values of residual volume and vapor ingestion time are given in 
Table I for the fourteen combinations of Weber number. Bond number, drain 
size, and initial filling. These data are also plotted as functions of the 
Weber and Bond number in Figures 21-23. Periods of draining induced waves 
observed to exist during the draining simulations are presented in Figure 2b. 

The program for simulating low gravity draining uses two new numerical 
methods that were developed specifically for it, namely, the method of variable 
trajectories and the method of moving triangular meshes. It has the property 
that most of the computation takes place at the free surface. Because cubic 
spline interpolation is used, few points are needed on the free surface. Al- 
though the free surface may more than double its length, no more than twenty- 
six points have been used. Computing times have varied from twenty minutes 
to an hour and three quarters, the latter .including nearly 8000 time steps, 
some rerunning, and sixteen oscillations. With additional work, the computing 
time can be reduced, by perhaps twenty percent. 
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With little modification the program could simulate the effect of throttling 
[7] and with some changes the effect of outlet baffles [8] on low gravity- 
draining at lower Weber numbers than can be achieved experimentally. An 
extension of the idea behind the method of truncated solutions can be used 
as a basis for three-dimensional draining problems with off-center drains. 

The program is also adaptable to computing simulations of axi symmetric, sur- 
face tension dominated reorientation problems. The methods can be general- 
ized to deal with the inviscid rotational approximation to the Navier-Stokes 
equations expressed in terms of stream function and vorticity, which is ap- 
plicable when the vorticity is produced outside the domain and merely con- 
served within. The methods can also be used for the full axisymmetric 
. Navier-Stokes equations. 

The low gravity draining program and its generalizations appear to have the 
same capabilities as the Marker and Cell family of programs [9] • The dis- 
tinction between the ranges of applicability lies in the philosophy behind 
each. The present program is based in classical continuum mechanics with 
classical mathematical techniques adapted for computer usage and achieves its 
results with extensive computation on few points. The Marker and Cell family 
is based upon doing a moderate amount of computation to imitate the physics 
in each cell and doing this for many cells and time steps. So long as the 
boundary holds together, or the domain breaks into large pieces, the family 
derived from the present program probably has an advantage; when the boundary 
decomposes irregularly, the advantage is to Marker and Cell. 

The numerical simulation program and the results obtained from it have interest 
beyond the field of low gravity draining. The demonstration that the combina- 
tion of the methods of variable trajectories, truncated direct solution, and 
moving triangular meshes can solve very difficult moving boundary problems 
over a wide range of parameters is a significant result to mathematicians 
interested in moving boundary problems. 
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PROBLEM FORMULATION 


In this section the mathematical formulation of the draining problem is 
presented. All the variables that appear are dimensionless, unless other- 
wise designated, and they are defined in the ’ Symbol List, Appendix A, where 
the relationships to the physical dimensioned variables are given. 


Initial Configuration 

Consider an axisymmetric vertical container consisting of a right circular 
cylinder with a hemispherical bottom that has an orifice of circular cross- 
section centered on the axis (see Figure 1(a)). Let the container be so 
oriented that the gravitational field is directed parallel to the axis. Prior 
to the initiation of draining at time t = 0, it is assumed that there is 
sufficient liquid in the container to cover the orifice and that the liquid 
is at rest. 


The equilibrium free surface of the liquid is found by solving the two-point 
boundary value problem 


with 


and 


r dF/dr 

(1 + (dF/dr) 2 ) 1 / 2 


1 _d_ 
r dr 

~ = 0 at r = 0 
dr 


- B(F - h ) = 2H , 0 < r < r 
v c' o’ w 


dF/dr 

(1 + (dF/dr) 2 ) 1 / 2 
at r = r w 


= r cos © - (1 - r) ' sin 0 


1/2 . 


, ( 1 ) 

( 2 ) 

( 3 ) 


Equation (l) is the requirement that the static liquid pressure at the free 
surface be in equilibrium with the constant ullage pressure there; (2) 
is the condition arising from the symmetry of the container; and (3) is the 
requirement that the free surface meet the container with contact angle 0 . 
(Computations are carried out in this study for the case where - 0 has the 
value 5 °)» Circular cylindrical coordinates are used with origin at the 



center of the orifice, and the Bond number B is chosen to be positive when 
gravity acts vertically downward (see Figure 1(a)). 

The volume of the liquid within the tank is 

r 

w 

- V = 2TT C r f(r) dr - V s (4) 

J o 

where, for an orifice of radius r , 

’ o' 


V 

s 


= Hf 
2 w 




is the volume between the hemispherical bottom of the tank and a right 
circular cylinder of radius r w with bottom touching the orifice (see Figure 
1(b)). By specifying, along with (l), (2), and (3), the initial volume V , 
the height F(r) of the initial equilibrium free surface is determined (and 
hence the values of h,, the height of the surface at its center, H , the 
mean curvature there, and r , the value of r where the free surface inter- 
sects the wall). 


Draining 

Dynamic Equations . At time t = 0, it is assumed that draining is initiated 
by pumping the fluid out in such a manner that the outflow velocity is directed 
vertically downward and is uniform over the orifice. This outflow velocity is 
attained instantaneously and remains constant until the time that the free 
surface first begins to enter the orifice. It is further assumed that the 
liquid is incompressible and inviscid and that vortexing and turbu- 
lence are absent so that the motion can be taken to be irrotational. In 
actual practice the effect of a small nonzero viscosity is evidenced by a 
layer of liquid that remains on the container wall during draining, but such 
a layer is generally thin enough for problems of interest here that it may 
be neglected. 

Under the above assumptions, one is led to the following initial boundary- 
value problem for axisymmetric motion from t = 0 until free surface 
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breakthrough, 


Velocity Potential . The sign of the velocity potential <p is chosen so 
that the velocity is given by v = ?cp . The velocity potential satisfies 
Laplace's equation. 



( 5 ) 


in the liquid. 


Boundary Conditions . The conditions on cp over the fixed portion of the 
boundary are that on the symmetry axis and container wall the normal velocity 
be zero 


I? 


= 0 


on r = 0 and on the wall 


and at the orifice the outward-directed normal velocity be 



( 6 ) 


dtp 1 

an = 2. 


for z = 0, 0<r<r 

o 


(7) 


On the instantaneous free surface, one has the kinematic condition, which 
can be expressed either in the form 


Sf _ clcp acp af 

at ” az ar ar 


dR _ a? dZ a 

dt a?' dt a dz 


(8a) 

(8b) 


and the Bernoulli equation, which can be expressed as either 

W 5t = 2H ‘ Bf - I M 2+ c OO (9a) 

or 
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wf2 = 2H-BZ + ||v<J)! 2 + C(t) 
at d — 


(9b) 



where 




is twice the mean curvature of the free surface, z = f(r,t). The first 
forms (8a) and (9a) are written in Eulerian (partial time derivative) form 
for points z = f(r,t) on a single valued free surface, whereas the second 
forms (8b) and (9b) are written in Lagrangian (total time derivative) form 
for points r = R(t), z = Z(t) on a free surface which may be multiply- 
valued (with 2H appropriately expressed). The first forms follow the motion 
of the free surface along fixed values of r, whereas the second follow the 
actual motion of points on the surface. 

The choice of the particular function of integration C(t) that, from 
theoretical considerations, is the most convenient to use in (9a) and (9b) 
is the one corresponding to cp being independent of time when the velocity 
at the free surface is uniform and directed vertically downward. For this 
case of stationary flow, we have that cp = const.- z near the free surface 
and f(r,t) = F(r) - t, so that setting d<p/8t = 0 in (9a) yields, with 
the use of (l), for C(t) the value, 

C(t) = -Bt + W/2 - (2 H q - Bh c ) (10) 


Finally, the end condition for (9a) and (9b), that the contact angle remain 
fixed during the motion at its equilibrium value ©, is 

= r cos 0 - (l - r^)^/^ sin 0 at r = r^ (ll), 


df/8r 

(1 + (df/dr) 2 ) 1 / 2 
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Initial Conditions . Equations (5) - (ll) are to be solved subject to the 
initial conditions at t = 0 that the free surface is the equilibrium shape 
2 as F(r) and that the velocity potential is the solution of (5), (6), (J), 
and the condition that the potential vanish on the free surface, 

tp = 0 on z = F(r) (12) 

That the free surface is initially an equipotential can be observed by 
examining the nature of the discontinuity in the Bernoulli equation (9) at 
t = 0 [lO]. Take c(t) equal to -2H + Bh for t < 0 and the value 
given by (10) for t > 0; then ( 9 ) holds for non-positive as well as for 
positive t with cp = const, for t < 0 . The right-hand side is piece- 
wise continuous with a step discontinuity at t = 0 corresponding to the 
instantaneous change of velocity. Thus the time derivative of cp at the 
free surface has only a step discontinuity at t = 0, and hence tp itself 
must be continuous there at t = 0 and equal to the same constant value it 
had for t < 0 . The continuity of <p on the free surface at t = 0 
corresponds to the condition that the external ullage pressure remain 
fixed there as draining is initiated [11 ]• 
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INC. 


NUMERICAL ANALYSIS 


This section begins with a summary of the numerical scheme used for simu- 
lating low gravity draining. After reading the first paragraphs in this 
part, the reader interested in the physical results of this study may Jump 
to the fourth part, which discusses the interpretation of data. The formal 
details of the new methods are covered in the second part. The third part 
describes some technical details of the program and includes, retrospec- 
tively, sketches of changes that would make the program even more success- 
ful. New problems that have arisen through the development of the program 
and through the consideration of its extension to other applications in 
low gravity fluid dynamics are also described. 

Numerical Scheme 

Reformulation of the Mathematical Problem . Vfith the free surface boundary 
conditions expressed in Lagrangian form, the mixed initial-boundary value 
problem becomes three sets of ordinary differential equations ( 8 b, 9 b) con- 
nected by the solution to a mixed boundary value problem ( 5 ), ( 6 ), and ( 7 ), 
for the velocity potential within the fluid, a domain with a moving boundary. 
Two of the sets of ordinary differential equations ( 8 b) describe the motion 
of particles lying on the meridian of the free surface; the third set ( 9 b) 
describes the change in velocity potential with time along the particle 
trajectories. Thus the mixed boundary value problem for the velocity poten- 
tial has given but changing boundary values on the free surface and, with 
constant drain rate, time-independent normal derivative conditions on the 
remainder of the boundary (cf. Figure l). This restatement of the mathe- 
matical problem was introduced by Moore and Perko [ 6 , p. 307 ] in studying 
reorientation in low gravity. 

The restatement implies that almost all the action is at the free surface in 
the following precise sense. The data at a point on the free surface needed 
to integrate the ordinary differential equations are the coordinates, the 
curvature and slope of the meridian, and the velocity potential and its 
derivatives, tangential and normal; of these only the last requires informa- 
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tion not on the free surface, and that, only in a neighborhood of it. In 
the nondimensionali zation the tank radius is taken as the characteristic 
length and the characteristic time is chosen so that the nondimensional 
outflow velocity integrates to 7T over the drain. Hence equal volumes 
leave the tank in equal times; and the drain radius remains as a geometric 
parameter controlling the shape of the flow through the drain. The dimen- 
sional outflow rate (and any variation in it) appear in the Weber number, 
which occurs only at the free surface in the ordinary differential equation 
for the change of velocity potential with time. (Were the drain rate to 
vary, the transformation back to dimensional time would vary too. ) The 
mixed boundary value problem has two functions; it transfers information 
along the free surface and it carries the residual geometric information 
left after nondimensionali zation up to the free surface. 

When the fluid is at rest, the initial shape is determined by balance of 
gravitational and surface tension forces, including the contact angle con- 
dition, and the free surface is initially an equipotential, which may be 
taken as zero. The computational realization of the restated problem has 
three parts: first, choosing a discrete set of points on the initial free 

surface as initial conditions from each of which to integrate numerically 
three ordinary differential equations; second, integrating the equations; 
and third, at each time step approximating the normal derivative of the 
velocity potential at the free surface by finding the potential along a line 
near and "parallel" to the free surface. A cubic spline interpolation gives 
a numerical representation of the free surface. 

The Method of Truncated Direct Solution . The statement, almost all the 
action is at the free surface, can be taken over to the discrete problem as 
most of the computation is near the free surface by taking advantage of the 
properties of the direct block factorization method [12, pp. 52-61 ] applied 
to the finite difference equations representing the mixed boundary value 
problem for the velocity potential. As yet only two restrictions on mesh 
and finite difference method are needed: first, that the mesh points can be 

arranged in lines (with those adjacent to the free surface roughly parallel 
to it), and second, that at each point the finite difference equation in- 
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volves only points on a few adjacent lines. With these restrictions, the 
totality of finite difference equations for the velocity potential can be 
arranged as a matrix equation in block form by lines with the blocks for 
the line adjacent to the free surface last. When a set of mesh lines 
adjacent to the free surface is changed, the last block rows (with differ- 
ence equations involving points in the changed rows) change too. However, 
because the boundary conditions, except on the free surface, are independent 
of time, the bulk of the matrix does not change; so the direct matrix does 
not change; so the direct matrix factorization down to the first changed 
block row can be saved, reused, and the process restarted there. Further- 
more, because only the potential on the line adjacent to the free surface is 
needed to approximate the normal derivative, the solution may, in general, 
be truncated after this line. If the vertical mesh spacing is wide enough 
near the free surface, the moving boundary will not come near the mesh line 
below it in a time step. Hence in the next time step only the last block 
row in the matrix equation need be recomputed and refactorized. 

Generation of Moving Triangular Meshes . As the simulation continues, it is 
necessary to revise the mesh lines to accommodate the changing shape of the 
free surface and the outflow of liquid. Change of shape is accommodated by 
revising the mesh progressively deeper and deeper into the tank; outflow is 
dealt with by removing a line from the mesh near the free surface. 

To facilitate approximating boundary conditions along curved boundaries, an 
irregular triangular mesh is used. Along the free surface two mesh triangles 
meet at the axial point, three meet at a general point, and there is one mesh 
triangle in the contact angle, in the interior of the liquid, six mesh tri- 
angles meet at a mesh point. With these restrictions, a triangular mesh 
filling the half-axial cross section of the liquid is topologically equival- 
ent to a polygonal domain composed of equilateral triangles, with the image 
of the free surface being a straight line (cf. Figure 2). The mesh lines 
within the liquid are the images of the lines parallel to the Image of the 
free surface in the polygonal domain. The generation of the triangular mesh, 
both initially and in revision, is accomplished by adapting a triangular 
mesh generator developed at Lawrence Radiation Laboratory, Livermore, 
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California p. 168] in which only the location of the boundary points 
need be prescribed. 

Finite Difference Equations . Each triangle in the mesh is subdivided into 
three quadrilaterals bounded by segments of the peipendicular bisectors of 
the sides when the triangle is acute and by segments of the medians when the 
triangle is obtuse. Thus each mesh point is surrounded by a polygonal area 
formed by as many quadrilaterals as there are triangles meeting at the mesh 
point. The union of these polygonal areas fills the half -axial cross section 
of the liquid. The finite difference approximation at each mesh point is 
obtained by integrating the Laplace equation over the mesh region generated by 
rotating the polygonal area surrounding the mesh point about the axis of 
symmetry. On one hand, this volume integral vanishes; on the other, it is 
equal to the integral of the normal derivative of the velocity potential over 
the bounding surfaces. The normal derivative is approxinated by the finite 
difference quotient between adjacent mesh points, unless it is given by a 
boundary condition. Thus each finite difference equation involves points on 
at most three lines. Ihe areas of , the bounding surfaces are sums of areas 
of the frustums of cones generated by the sides of mesh regions. This pro- 
cedure gives the finite difference equations the weights appropriate to the 
axial symmetry of the problem while retaining the advantage of calculating 
mainly in two dimensions. Additional details are given throughout [l] and 
in [ 13 , p. 151; 14, pp. 181-194] . 

Computation of the Derivatives . For calculating the curvature and slope of 
the free surface meridian, its two coordinates are interpolated by cubic 
splines in chord length along the polygonal free surface. These interpolants 
have end values implied by the fixed contact angle and the symmetry at the 
center, pass through the computed values, and have continuous first and 
second derivatives [ 15 , p. 51] • This permits accurate calculation .of the 
curvature from few mesh points on the free surface and accurate rotation of 
the normal and tangential derivatives into the derivatives' i n the coordinate 
directions. The tangential derivative of the velocity potential is deter- 
mined from a similar cubic spline interpolation; the end conditions are the 
zero derivative implied by symmetry at the center and a cubic extrapolation 
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of the derivative to the wall. At each mesh point on the free surface, the 
integral of the normal derivative over the bounding surfaces of the corres- 
ponding mesh region must vanish. Solving the resulting finite difference 
approximation for the unknown normal derivative gives the latter as a linear 
combination of five values of the velocity potential, three on the free sur- 
face and two on the mesh line below. Because the contact point must be con- 
strained to move along the wall, its velocity is taken as the derivative of 
the velocity potential along the chord of the wall. 

Normalization of the Potential and Integration . For some cases, particu- 
larly those with low Weber number, the potential may become very flat near 
the contact point, lb minimize the loss of significance in forming differ- 
ences, the arbitrary function of time, C(t), in the Bernoulli equation, (9a) 
or (9b), is chosen so that throughout the computation the velocity potential 
is zero at the contact point. Because the derivative of the velocity poten- 
tial must vanish there, C(t) is the negative of the other terms in the dif- 
ferential equation for the change of velocity potential at the contact point. 
When the contact point lies on the hemisphere, it is returned to the hemi- 
sphere by projection from its center at the end of each time step. The 
trapezoidal rule is used for the numerical integration of the ordinary dif- 
ferential equations. Time steps are varied slowly, in geometric progressions 
with ratios near one; the user prescribes the interval for which each ratio 
holds. 

The Auxiliary Potential. The scheme for simulation so far described still has 
an unresolved numerical and programming problem. Relative to the mesh that 
is adequate at the free surface and throughout the tank, a fine mesh is needed 
to compute the flow into the drain. To get reasonable residuals, the free 
surface must approach, if not enter, this region of fine mesh. Thus it ap- 
pears that mesh points should be periodically added to the free surface with- 
out introducing unwanted oscillations thereon. To avoid this problem, which 
is probably lengthy but not insoluable, an auxiliary potential is introduced 
to satisfy the inhomogeneous boundary condition on the velocity potential at 
the drain. The auxiliary potential is that of the same liquid draining at 
the same rate from a circumscribing, flat-bottomed, semi-infinite, circular, 
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cylindrical tank through a circular, axial drain of the same size (cf. Figure 
l). When the potential is represented as the sum of the auxiliary potential 
and a residual potential, the latter satisfies. a mixed boundary value prob- 
lem which has boundary conditions identical in form to those for the full 
velocity potential except that the normal derivative vanishes in the drain 
and is a nonzero function on the circular part of the tank meridian. This 
function represents the negative of the flow into and out of the corner of 
the cylinder cut off by the tank meridian; it is continuous and slowly vary- 
ing and matches the zero normal derivatives it meets at its end points. 

Hence it places no undue restriction on the mesh. Introducing the auxiliary 
potential transfers the major part of the flow geometry at the drain to the 
free surface. In spite of its drawbacks, extra computing time and extreme 
difficulty in evaluating near the drain, the introduction of the auxiliary 
potential was the key that made possible the realization of the numerical 
scheme Just described in the previous study. The auxiliary potential has 
been carried over to the present study. 

The Method of Variable Trajectories . Adding a small auxiliary velocity to 
the tangential velocity on the free surface meridian and leaving the normal 
velocity unchanged result in a small first-order displacement of the compu- 
tational point along the free surface at the new time step. Repeatedly doing 
so gives variable integration trajectories. Because the axis and the wall 
are fixed trajectories, the auxiliary velocity must be restricted to vanish 
on these lines. Choosing the auxiliary velocity identically zero gives the 
Lagrangian formulation of the problem; choosing it so that the trajectories 
lie along the radial coordinate lines results in the Eulerian formulation. 
Because the free surface stretches and contracts nonuniformly and contains a 
fixed number of mesh points, the rule used for choosing the auxiliary velocity 
is one that maintains the ratio between adjacent mesh intervals nearly con- 
stant. The differential equation for the change in velocity is simply the 
total derivative with respect to time along the trajectory, which can be com- 
puted from the Bernoulli equation (9 a )> the sum of the tangential and the 
auxiliary velocity, and the normal velocity. The introduction of variable 
trajectories prevents the growth of large mesh spacings during the development 
of a wall sheet and the movement of mesh points too close to the axis. 
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The Method of Moving Triangular Meshes . The distance between two mesh lines 
is directly available only at the centerline and the wall; however, because 
it Is the altitude of mesh triangles touching the two mesh lines, it can be 
inferred from the area of the triangles, which is always available in the 
calculation. Mesh triangles touching the free surface are divided into four 
groups of contiguous triangles, upper and lower bounds for the areas within 
each group are calculated at each mesh revision, and these bounds are ex- 
tended to test bounds by multiplying by fractions near one, which may vary 
from group to group. At each subsequent time step the areas of triangles 
touching the free surface are tested against group bounds. When the area of 
a triangle falls outside the test bounds, rezoning is initiated. Thus re- 
zoning takes place when the distance between the free surface- and the adjacent 
mesh line becomes too small or too large. When wave motion is strong, rezon- 
ing may be initiated from any group; when a wall sheet is developing, gener- 
ally it is from the center line group. 

The first rezoning is a minor one involving three rows of mesh triangles ad- 
jacent to the free surface. At the first minor rezoning, test bounds similar 
to those used at the free surface are set up for the lowest row. When minor 
rezoning is again required, the areas in the lowest row are checked against 
the test bounds. Failure indicates that some part of the top three rows has 
been compressed (or dilated) to the point where readjustment of the mesh over 
more rows is necessary. A major rezoning over five or six mesh rows from the 
top frequently suffices. Again the first time a major rezoning occurs, test 
bounds are set up for the lowest row of triangles and the sequence of minor 
rezoning tests is restarted. At the next major rezoning, the areas of the 
lowest row of triangles are tested against the test bounds. ' Failure requires 
a general rezoning deeper into the tank. The test for determining when re- 
zoning is necessary is carried over from the previous study; the analogous 
tests for determining whether the bottomline of a mesh revision are satis- 
factory are new. Their failure is a prediction that further rezoning will 
soon be required and an indication that an immediate -rezoning will save time. 

The mesh generator requires that the mesh points on the boundary of the region 
to be rezoned be prescribed. Because the mesh points on the free surface and 
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lowest mesh line are fixed, the rezoning input reduces to specifying the mesh 
points on the axis and wall within the region to be rezoned, lb insure that 
the mesh varies smoothly, particularly near the free surface, these points 
are distributed in geometric progressions whose ratios are determined by the 
number of intervals, the initial interval, and the length of the rezoning 
interval. At the wall the first interval is set by making the triangle at 
the contact point isosceles; at the axis it is a fraction of the adjacent mesh 
spacing at the free surface. 

Before anything else is done, at any stage of rezoning, these ratios are com- 
puted and tested against input test bounds which may differ at the wall and 
the axis and usually become more stringent as rezoning progresses from minor 
to major to general rezoning. Failure to pass the test bounds is cause to 
progress to the next stage. Passing the test bounds for minor and major re- 
zonings invokes the area tests described above. In the general case two more 
mesh rows are rezoned after a failure. After successive failures to satisfy 
the test bounds, a mesh line is dropped between the free surface and the 
current fixed mesh line. When the general rezoning reaches the bottom of the 
tank, the final mesh row is adjusted, if necessary, to meet the ratio require- 
ments along the wall, all switches are reset to restart the rezoning process 
from the beginning, and the entire mesh is rezoned. Figure 3 shows selections 
from the sequence of meshes generated for an oscillatory problem. 

The rezoning procedure just described is satisfactory so long as a marked wall 
sheet does not develop. When the fourth mesh point on the free surface, 
counting the contact point as the first, moves toward the wall across a test 
bound, .958 has been used, programming is Invoked which removes all but the 
top four triangles in the right-hand diagonal column from the logical diagram 
and the corresponding partial column from the physical mesh (cf. Figure 2). 
When the next point on the free surface moves across the test bound toward the 
wall, the next diagonal column is dropped, except for the top four triangles. 
The dropping continues so long as mesh points on the' free surface cross the 
test bound. Of course, there is also programming to restore these columns to 
the mesh as mesh points cross the test bound moving away from the wall. Thus 
the logical diagram for a developing wall sheet is the union of two parallelo- 
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grams having a common left boundary with the wider Just below the free sur- 
face containing two rows of triangles. Mesh points on the right and lower 
boundaries of the logical diagram map into the tank meridian in the physical 
mesh. The radial coordinate oi tne points which lie on the exterior boundary 
of the top parallelogram is one; their vertical coordinates are the height 
at the free surface of the corresponding points on the free surface counting 
from the contact point. The last mesh interval so determined plays the role 
of the side of the isosceles triangle as the first mesh interval in generating 
subsequent rezoning boundary points along the wall. 

The number of mesh intervals remaining on the right edge of the logical dia- 
gram is known. The length of the interval to be filled with a geometric 
sequence of mesh intervals is determined by assigning the physical mesh point 
corresponding to the lower right corner of the old logical diagram to the 
same corner in the new one. Along the arc corresponding to the bottom of the 
logical diagram, less the intervals mapping into the drain, mesh points are 
distributed uniformly. Along the axis mesh points are distributed uniformly 
between free surface and drain to complete the specification of boundary 
points necessary for a complete rezoning. 

As the wall sheet develops, the center falls and mesh lines are dropped. Hence 
with few mesh lines, the ratio of the geometric sequence along the wall in- 
creases and may exceed the upper bound prescribed for a general rezoning. 

When this occurs, another mesh interval is added on the wall by dropping part 
of another diagonal mesh column. The lower part of the right boundary of the 
logical diagram is moved one mesh interval to the left, leaving six mesh tri- 
angles unchanged. The logical diagram is now a union of three parallelograms 
growing wider toward the free surface. The heights of mesh points on the wall 
that correspond to the exterior boundaries of the top two parallelograms are 
obtained by projecting the height of the free surface at the corresponding 
point counting from the contact point along both edges. Mesh intervals between 
boundary points corresponding to the lower exterior edges of the top parallelo- 
grams are extra mesh intervals on the wall, in that they do not correspond to 
mesh intervals at the axis. Ihe length of the last extra mesh interval is 
again used as first interval in the geometric sequence below. Complete re- 
zoning as described in the previous paragraph takes place. 
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'Hie program permits adding extra intervals of either kind as needed. Figure 
2 shows a physical mesh and the corresponding logical diagram with seven 
extra mesh intervals, four generated by motion of mesh points toward the wall 
and three required to keep the mesh intervals along the wall from becoming 
too large. 

Extra mesh intervals generated by the motion of points on the free surface 
are retained so long as they lie completely on the vertical tank wall and the 
corresponding point on the free surface lies to the right of the test bound. 
When the latter condition fails, for one or more points, the corresponding 
intervals are removed and the outer boundary of the logical diagram below the 
second row is moved to the right. When any extra interval begins to escape 
from the tank wall, it is removed and the corresponding mesh column is com- 
pletely restored to the logical di agism and the physical mesh. On the occas- 
ion of the addition or deletion of an extra mesh interval along the wall, the 
entire mesh is revised. 

Computational Normalizations and Checks . In the program the definitions of 
the volume, the rate of flow out the drain, and the streamfunction are divided 
by n . With this modification the streamfunction varies between 0 at the wall 
and -1 at the axis (or 0 at the axis and 1 at the wall, If the other direction 
of Integration is chosen). Thus the streamlines can be naturally labeled so 
that a given streamline is the trace in the half -axial cross section of the 
curved cylinder enclosing the given fraction of the flow out the drain. This 
labeling facilitates interpreting the plots showing trajectories of stream- 
line ends at the free surface. Another consequence of the renormalization is 
that the modified rate of flow out the drain is unity and the total modified 
outflow in a given time is equal to the time. Thus the principle that the 
sum of the volume in the tank and the volume that has left the tank should 
equal the original volume translates to the statement that the computed volume 
plus the time should equal the initial computed volume. Volumes are computed 
from the cubic spline representation of the meridian of the free surface by 
Lobatto quadrature of order six. 

The value of the streamfunction at the axis is printed at each time step. A 
drift away or a jump from a very good approximation to -1 indicates diffi- 
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culty with the mesh or with the computation of the auxiliary potential. As 
an over-all check of the accuracy of the computation, the sum of the time 
and the volume is printed at each time step. 

Small Amplitude Waves. Because there is no requirement that the free surface 
be single valued, the method of variable trajectories shares with the 
lAgrangian formulation the ability to detect waves of small amplitude, par-, 
ticularly along a wall sheet, provided that the mesh spacing is fine enough. 
(This does occur in Figure 13(c) and (d).) However, using a small mesh not 
only materially increases computing time but also may lead to computational 
instabilities that defeat the goal of determining the over-all behavior of 
low gravity draining. 

New Methods 

The Method of Variable Trajectories . Here the formal details of the method 
are given (its use as part of the numerical scheme is described above under 
the same heading). Throughout the computation, a local coordinate system is 
used at the free surface. The exterior normal and the tangent vector di- 
rected toward the axis are the local unit vectors. The method of variable 
trajectories moves points on the free surface by adding to the tangential 
velocity a small incremental velocity, u, also directed along the free sur- 
face. For the current problem the function u is required to vanish at the 
axis and the wall. The equations for the evolution of the free surface, the 
analogs of the kinematic equations (8b), take the simple form 

N , = cp and S , = cp + u (13) 

t n t s 

with N(t) and S(t) representing position in local coordinates. The 
corresponding equation for the evolution of the potential is the total deriva- 
tive of the potential along the variable trajectory, 

wfj = 2H - BZ + - | V cp | 2 + u cp„ + C(t) (1*0 

at d 5 

where 2H = (Z R - Z R ) + Z /R (l5) 

ss s s ss s' 

is the mean curvature of the free surface [l6, pp. 7> 11 ] . The transforma- 
tion of (13)' from the local coordinates to the global cylindrical coordinates 
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is accomplished by a rotation through the angle, - arctan (Z /R ), followed by 
a reflection in the r-coordinate line through the point, which occurs because 
the tangent vector is directed away from the axis in the global system. The 
resulting computational equations for the evolution of the free surface 

R, = -cp Z - (cp + u) R 
t n s s 

Z. = <p R - (cp + u) Z ^ 

t ns ' s s 

differ from their Lagrangian counterparts only in the incremental velocity. 
Aside from the choice of u, only one multiplication and two additions are 
added by (l 4 ) and (l 6 ) at each point. These equations are used in both the 
predictor and corrector steps of the modified Euler method. 

The point obtained by applying only the second equation in (13) over a 
finite time interval does not lie on che free surface. This observation sug- 
gests that a second order auxiliary velocity, based upon the curvature of the 
free surface (the term in parenthesis on the right of (15)), could be added 
to the normal velocity to give a second order correction to the trajectories. 
Although the information needed is available, the first order modification 
proved adequate for the present applications. 

The goal of the choice of u is to maintain the ratio of adjacent mesh inter- 
vals in the neighborhood of a fixed input value. Because the number of 
intervals in the geometric sequence is given, the size of the initial interval 
varies in response to the stretching and contracting of the free surface. 

From this information, a geometric sequence of points with the given ratio 
between adjacent intervals is determined. The difference in chord length 
along the free surface between the actual mesh point and the corresponding 
member of the desired geometric sequence is reduced to the velocity, u, by 
dividing by five times the current time step. The multiplier five, which 
is arbitrary, and the further restriction that the magnitude of u not exceed 
that of the tangential velocity are devices, adopted in- lieu of a more com- 
plicated rule for determining u, to insure that the trajectories near the 
wall remain within the tank. These modifications are appropriate because 
continued progress toward the goal is what is needed. In comparison with 
its effect, the additional work in determining this u is modest; five addi- 
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tions and four multiplications per point plus the evaluation of a logarithm 
and an exponential and some overhead are needed. 

Experience gained in this study suggests that replacing the auxiliary 
potential by a finite difference approximation throughout the tank, with 
many mesh lines running into the drain, may be feasible. In this case mesh 
points would have to be added to the free surface as it approached the drain. 
The method of variable trajectories, with an appropriate choice of u, could 
be used to spread mesh points apart In anticipation of adding a point or 
points in regions where unwanted instabilities would not develop. 

The Method of Moving Triangular Meshes . The application of the method as 
part of the numerical scheme is described above under the same heading; here 
further details and some of the underlying principles are discussed. A non- 
uniform triangular mesh in which six mesh triangles meet in an interior point 
can be deformed into a regular triangular mesh with a polygonal boundary, 
the corresponding logical diagram. Frequently, the discussion of the struc- 
ture, generation, movement, and change of a nonuniform triangular mesh be- 
comes simpler in terms of its logical diagram. 

For programming purposes, the polygon in the logical domain should be the 
union of simple geometric figures composed of mesh triangles, such as tri- 
angles, parallelograms, and trapezoids. And it should have a connected 
Interior. The mesh triangles within the polygon between two adjacent mesh 
lines form a mesh row. The programming is simpler when all the mesh rows 
in at least one of the three directions are simply connected. At the cost 
of extra space and some programming, multiply connected domains can be 
handled with existing coding. For numerical reasons, the mesh generator 
mapping the logical into the physical mesh should require minimal input, vary 
the interior mesh smoothly when the mesh on the boundary is smooth, and pro- 
duce interior triangles that are as nearly regular as possible. The last re- 
quirement is to insure that no direction is slighted in the flow of informa- 
tion within the mesh. Obvious violations of these requirements occur in the 
present study in the neighborhood of the wall sheet and above the drain. 
Experience has shown that when the free surface is close to the drain, 
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too little mesh in the vertical direction either produces an oscillation or 
enlarges a real one. Further numerical experiment with a finer vertical 
mesh is needed to determine the nature of this oscillation. 

In the physical mesh, the free surface is a distinguished curve along which 
it is desired to keep the mesh structure as uniform as possible. Thus the 
same number of mesh triangles should meet at each mesh point. Hence the 
image of this curve in the logical diagram should be a straight line. The 
method of variable trajectories is designed to enforce a measure of regu- 
larity with a fixed number of mesh points, even when the free surface is 
stretching and contracting. The succession of minor, major, and general 
rezonings. are devices to restore uniformity near the free surface with a 
minimum of work. Stability of the normal derivative on the free surface at 
the axis is a continuing problem. Two devices have minimized this problem; 
first, the number of rows in a minor rezoning was increased from two to 
three, and second, the axial mesh spacing was made uniform. In major and 
general rezonings, the aim is to relax the compression and dilation near the 
free surface progressively further into the mesh. Uniformity along the 
axis near the free surface is achieved by requiring that the ratio between 
successive axial mesh intervals be near one for a major rezoning and even 
nearer for a general rezoning. 

In changing the logical diagram to accommodate the development of the wall 
sheet, the changes were made as far away from the moving boundary as appeared 
practical. Experience shows that reliable results can be obtained with quite 
irregular meshes far from the free surface. More mesh points on the free 
surface are desirable near a base of the thinning wall sheet. In this area 
there may be an abrupt transition very near the free surface between long thin 
triangles in the sheet and the fat triangles in the body of the liquid (cf. 
Figure 16(f)). It would be useful to know under whst conditions the cubic 
spline interpolation could be used to readjust and add mesh points along the 
free surface without triggering computational instabilities. 

The method of moving triangular meshes appears to be applicable to a number 
of moving boundary problems. Among those in low gravity fluid dynamics are 
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the investigation of the effect of baffles in the bottom of the tank and 
various aspects of reorientation including the generation and dissipation 
of eddies and geysers. 

Details of the Program 

Only the aspects of the program relevant to the understanding of the present 
simulations are discussed here. Material not covered here or in the earlier 
parts of this report is unchanged from [l 3 . After a description of the 
present organization, a few specific details are given in the order that they 
occur in the first pass through the program. 

Organization . The program is a special purpose code for numerical simulation 
of low gravity draining; the experimental nature of the problem requires a 
highly modular design for ease in modification. Thus the program can be 
easily adapted for related problems in low gravity fluid dynamics. In com- 
plexity it is comparable with members of the Marker and Cell family and with 
large structural dynamics codes, although it uses only 56,000 word6 of core 
storage. It generates six files on auxiliary disk storage, one containing 
restart points, two for intermediate output, and three for storage of output 
for future processing, plotting, summary listing, etc. It is designed for a 
time-sharing computer system which, in the main, charges for the number of 
operations performed and for the product of the amount of core memory used 
times the time it is occupied. Therefore, overlays are used insofar as pos- 
sible; and subcollections are used, with the dimensions of arrays and buffers 
set only in the final collection of the absolute program. 

Generation of the Free Surface. For Bond number 0, an arc of a circle making 
the given contact angle can be used as an initial free surface meridian for 
experiments with contact angles less than 5°* To go to contact angle zero 
would probably require an analytic expression as a patch at the contact point. 
For higher Bond numbers, the program for integrating (l), (2), and (3) should 
be revised for smaller contact angles or for fill levels intersecting the 
wall within the hemisphere. An interesting and perhaps more instructive 
alternate would be to by pass the generation by simulating a full experimental 
run including the reorientation phase. 
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Initial and Restart Input . Two forms of restart are available. One con- 
tinues with the mesh from the previous part of the run; the other permits 
the user to improve the mesh. In both cases all control parameters may be 
modified. 


Matrix Generation . Were a reorientation phase to precede the draining phase, 
the change would come here. At the instant the drain opened, there would 
be a switch from the zero normal derivative conditions along the entire tank 
wall and in the drain to the nonzero conditions along the hemisphere represent- 
ing the residual influence of the auxiliary potential. The impulsive velocity 
resulting from the opening of the drain would be added to the existing velocity. 

Matrix Factorization . The chief reason for using few mesh points is the rapid 
increase with the number of points of computer time used in matrix factoriza- 
tion together with the solution of the lover triangular system of equations. 
These two subroutines generally use more computer time than any other part 
of the program. These routines deal with block- tridiagonal matrices having 
diagonal blocks of varying size. Inner products are computed In double pre- 
cision and truncated to single precision to strike a balance between accuracy 
and storage. 

Auxiliary Potential . For purposes of numerical computation, it is convenient 
to introduce the residual velocity potential tp^, 

<P 2 = <P - (17) 

where cp^ is the auxiliary potential corresponding to the draining problem in 
a full infinitely- tall cylinder having a flat, rather than hemispherical, 
bottom. This, in essence, transfers the nonhoraogeneous boundary condition 
on <p at the orifice to one on cp^ at the free surface, thus obviating the need 
for an accurate -numerical solution near the drain. The potential <p^ can be 
expressed as the infinite series 
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( 18 ) 


where and J, are Bessel functions of the first kind and J.. is the m-th 
u 1 1, in 

zero of the latter (taking Q = 0). For computational convenience the 
additive constant c is initially chosen so that cp^ = 0 at the intersection 
of the equilibrium free surface and the container wall. Substitution of 


V 


9^ for V into the draining equations yields the desired problem with 


which does not vary rapidly near the orifice, as the unknown. (See 
Figure 1(b)). 


Whenever the mesh is completely rezoned, the additive constant c is redefined 
so that the auxiliary potential vanishes at the contact point. High in the 
tank, where z is large, the convergence of the series (l8) is enforced by 
the rapid decrease of the exponential factor. As z decrease^ the converg- 
ence becomes slower and slower; and the computation of the auxiliary potential 
becomes a dominant factor in the computing time for low Weber numbers. At 
the high Weber numbers used in the previous study, the draining time is so 
short that this fact was obscured. An improved method of calculating the 
auxiliary potential is needed. A more rewarding alternate would be to elim- 
inate the auxiliary potential altogether. For details of the present calcu- 
lation see [l, pp. 19 - 21 ] . 


Treatment of the Contact Point . Because the contact point must be constrained 
to move along the wall, its velocity is taken to be the difference quotient 
of the potential at the contact point and the mesh point Just below on the 
wall. This choice is consistent with the calculation of Hie normal derivative 
at other points on the free surface. It also enforces on the vertical wall 
the boundary condition that the normal derivative of the potential vanish and 
provides a reasonable approximation to this condition on the curved wall. The 
value of 2H used in (l4) is computed at the contact point Just as at any 
other point. The condition that the contact angle remain fixed appears in 
the calculation as an end condition for the cubic spline interpolation. This 
is the only place where the contact angle condition occurs in the body of the 
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program. The successful computation with a very thin wall sheet suggests 
that a smaller contact angle may he feasible. The extrapolated value of 
the tangential derivative at the contact point is used in the calculation 
only as an end condition for the cubic spline interpolation. All 
values used at the contact point can be viewed as interior limits along 
the free surface or wall. 

Integration Control . To achieve stability in a vide range of situations, 
the modified Euler method is used for the integration of the differential 
equations for the location of the free surface and the potential thereon. 

The dependent variables are carried in double precision and the derivatives 
in single precision. Time steps are taken in geometric series with the 
ratios and times of ratio change under user control. What is needed here, 
as in many other problems. Is a simple and effective method for automatic 
time step control for low order methods where the first neglected term may 
not adequately represent the error. 

Mesh Control . One of the achievements of this study is to have demonstrated 
that the method of moving triangular meshes can be made to work effectively. 
The controls have already been described in the two preceding parts of this 
section with the method. 

Use of Plots and Data 

Plots . Streamline plots presented in this report are truncated after the 
data for the first ten mesh rows below the free surface has been plotted. As 
the previous study showed that the streamlines in the middle and bottom of 
the tank remain fixed [l, pp. 6 3 , 67 , 95 1 > continuing the plots beyond the 
truncation would waste computer time in redetermining information available 
on plots produced later in the run. Because the integral of the nondimen- 
sional outflow velocity over the drain is normalized in the computation to 
the value +1, the normalized streamfunction varies between 0 and -1 when the 
direction of integration is from wall to axis. With this normalization the 
streamlines can be labeled so that -.5 streamline is the trace (in the half- 
axial cross section of the tank) of the curved circular cylinder that instan- 
taneously separates the flow through the tank and out the drain into two equal 
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PAGES MISSING FROM AVAILABLE VERSION 


using the height at the centerline or wall as a common entry in the two 
graphs. “Hie columns on the left, headed STEP and TIME, give the time-step 
number and the time for each free surface included in the plot. Alternate 
entries in the columns are offset as an aid in counting. Within a simula- 
tion run, the time step may vary and the number of time steps between free 
surface output may either change or be offset at a restart. The free sur- 
faces presented here are an automatic selection from the original data chosen 
to give readability on a single page; consequently, intervals in step number 
and time in the lists may be unequal. 

In the labels on computer plots, the symbol Cl means Case lj B and W have 
their usual meaning; VOL is the nondimensional volume divided by7r; and DR 
stands for the nondimensional drain radius. Labels similar to those on the 
free surface plots appear at the top of the mesh plots with the corresponding 
values at the bottom. Small circles at the center of a triangle indicate 
that it is obtuse (an obtuse triangle at the intersection of the axis with the 
drain is a digital artifact produced by roundoff). 

Termination . In Table I, Summary of Computed Data, the values are from the 
printed computer listing of the last completed time step. For later simula- 
tions, the run was stopped when the free surface fell below .06 at the axis; 
some of the earlier simulations were terminated higher in the tank when it 
was clear that a restart would run for only a few more time steps and produce 
no significant change in the values. The data files (now on magnetic tape) 
may contain a free surface position from the predictor half of the next time 
step. • This predicted information appears as the last point of the time 
histories of the centerline and wall heists. An *on these plots marks the 
time when the predicted free surface jumped from ab'ove .0 6 to below the ori- 
fice of the drain. If the predicted free surface survives the automatic 
selection process, it gives a linear interpolation of 1 he mesh points on the 
free surface. The streamline end trajectories are computed and plotted for 
those time steps at which a complete free surface is not output. As vapor 
ingestion develops, the trajectories become separate dots near the axis. 

The dots are computed from the cubic spline interpolation of the free sur- 
face and therefore should be connected by a smooth curve perpendicular to 
the axis. 

32 


LOCKHEED MISSILES & SPACE COMPANY. INC. 


Determl nation of the Periods . When the computational phase of the study- 
terminated, the data available for the determination of the periods was a 
set of computer printouts listing h c , v fi , h^, and v^, the heights and vel- 
ocities at the centerline and wall. If the centerline were falling at a 
constant rate and had a purely periodic component U(t) with period A, then 
it would have the form h^ = a+bt+U for constants a and b; and the velocity 
would take the form v^ = b+U , a periodic function displaced by an unknown 
constant b. With this hypothesis, the maxima and minima of should, 
separately, occur with period In fact, the strongest maxima and 

minima of v are separated by intervals that tend to cluster about a common 
c 

value. Some of the zeros of v , extrema of h c , are also separated by values 
that belong to the cluster. A similar analysis applied to the functions h 

w 

and should yield an overlapping cluster; and this is the case. The values 
of the periods reported in Table II are averages of estimates determined by 
this procedure. 

The data used in this analysis may be complicated by the special nature of 
the centerline and wall and certainly are confused by nonlinear wave phenomena 
and by effects of the hemispherical bottom. Although more elaborate statis- 
tical and analytic schemes could be used to determine the periods, their ap- 
plication seems unwarranted until the nonlinear feedback and interference 
phenomena are better understood. In retrospect it would have been desirable . 
to integrate the Eulerian form of the equations, (8a) and (9a), along several 
fixed radii. This would increase confidence in the determination of the 
periods and provide additional data on wave motion. For studying the latter, 
detailed information on the free surface curvature and the potential and 
kinetic energy terms of the Bernoulli equation along the free surface would 
be desirable for selected portions of runs where phenomena of interest have 
been observed. Plotting trajectories of points on the free surface where 
these functions have constant value on a graph showing successive positions of 
the free surface would be a simple way to organize and to begin the analysis 
of this data. 
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Mondimensionalization and Normalization. The nondimensionalizing relations 


used in this study are 

X = aX, t ~ t7Ta 3 /C, and <p = (pQ/rra. 


(19) 


where the quantities topped by a bar are dimensional, a is the characteristic 
length (tank radius), 0 is the volumetric drain rate, and X is a generic 
spatial variable. Here the unit of nondimensional time is that time needed 
to drain a right cylinder of radius a and height a at a given drain rate. 

This time scale may be inappropriate when draining excites axisymmetric oscil- 
lations of the free surface. The nondimensionalizing relations used in 


previous sloshing work 


17, p. il] are an alternative. They are 


X = aX, t = t (l+B) o /( p a^) ” 2 , and <P = cp (l+B) aa/p 


( 20 ) 


where t and (p are nondimensional sloshing time and potential, B is the 
Bond number, and p and a are the dimensional liquid density and surface 
tension. The relations connecting the two nondimensionalizations take the 
form 

' t = t*K and (p = <p */ K with K = (w/(l+B)) 2 (2l) 


and show that the two nondimensionalizations are equivalent when K is one. 

For values of the ratio K less than one, the principal concern in developing 
the program for this study, the draining potential is the larger of the two; 
and the draining time is the smaller. One justification for continuing the 
draining nondimensional! zation from the previous study is that the greater 
potential range for K less than one facilitates accurate determination of 
potential differences, a continuing source of difficulty in the program. For 
values of K greater than one, the draining nondimensionali zation was known to 
be appropriate. However, the sloshing nondimensionalization (20) and its re- 
lation to the draining nondimensionalization (2l) play an important theoretical 
role in discussing the results. 

Accuracy . The principal check on overall accuracy is the conservation of 
volume. The average percentage change in volume from start to finish for the 
last eight cases run is .02h% with a spread from .001 to .065. Cases 1 and 2 
(see Table I for parameters) have a change of less than . 276$, which will be 


3 ^ 


LOCKHEED MISSILES & SPACE COMPANY. INC. 


discussed in the next section. Time does not permit extracting this data 
from the early runs, Cases 5> 29, 32, and 35} but no cause for concern was 
observed in earlier analysis. The value of the stream function on the free 
surface at the axis approximates -1 to four or more digits, except that two 
or three digits may be acceptable for very long running problems very low in 
the tank. 

Irregularities in the streamlines, over and above those due to digitalization 
in plotting, are visible in the comers at the free surface and just above 
the drain. Those at the free surface are the result of poor determination of 
differences of the potential in regions where it has a small gradient. Those 
near the drain result from the difficulty in evaluating the auxiliary poten- 
tial low in the tank and are sometimes compounded by too coarse a mesh. 
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DISCUSSION OF RESULTS 


Simulations of low gravity draining have been computed covering a wide range 
of parameters: two Bond numbers, 0 and 5; three and one half decades in 

Weber number, .01 through 50; two initial volumes, 5 7t/3 and 8 7r/ 3 (flat 
filling in normal gravity to two and three tank radii); and two drain radii, 
l/30 and l/lO of the tank radius. Fourteen simulations have been computed 
and the data is presented graphically in Figures 4 through 17. 

Each figure begins with two facing plots which provide a two-dimensional 
representation of the three-dimensional time history of the simulation. 
Subsequent pages may give expanded scale plots or streamline patterns illus- 
trating significant details of the simulation. These data figures are ar- 
ranged, for ease in cross reference, with those for Bond number 0 preceding 
those for Bond number 5 and within each group in order of increasing Weber 
number with the higher fill levels last. Each simulation is identified by 
a case number; and its parameters are identified by a symbol (B, W, r Q , V q ); 
so the expression Case 1 (0, .01, l/lO, 5 ^/3) means that Bond number 0, 

Weber number .01, drain radius l/lO, and initial volume 5 n/3 are parameters 
for Case 1. Table I contains a list of the case numbers and the correspond- 
ing parameters for ready reference. 

Because the simulations are solutions of an initial boundary value problem 
that is well posed both numerically and mathematically, their physical 
features should vary continuously with the parameters. This principle is 
the basis for the analysis in this section; however, with fourteen data 
points widely spread in a four-parameter input matrix, identifying continuous 
variation is sometimes difficult. This study pinpoints some open questions. 

The analysis of low gravity draining phenomena given in the first part of 
this section is based upon the totality of simulation data; however, it is 
presented when possible in terms of the graphical data in Figures b through 
17. Next, two simulations are compared with two experiments performed at 
NASA- Lewis Research Center. The agreement is excellent; on the basis of the 
preceding analysis, many common details are revealed. 
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An analysis of the numerical data, presented in Table I, follows. Ihe liquid 
residual volume, V^, is expressed as a fraction of a hemispherical volume; 
the initial volume, V q , is expressed throughout this report as nondimensional 
geometric volume. Liquid residual volume is discussed as a function of the 
problem parameters and the correlation variable W/(l+B). Graphical summaries 
of the data from this study and its predecessor [l] are presented in form 
suitable for engineering use. The final part of this section concerns wave 
motion. It is shown that the period of the axisymmetric oscillations excited 
by draining, when normalized to the sloshing time scale, is a constant charac- 
teristic of Bond number. 

Analysis of the Simulations 

As background for analyzing the data presented in Figures k through IT, the 
wall rise and initial transient phenomena are discussed. In the analysis, 
cases are grouped according to common physical phenomena or a common parameter 
change. Because there are four continuous parameters to vary, any single 
classification scheme is unlikely to be satisfactory. However, most schemes 
would include a capillarity dominated regime characterized by low Weber 
number, low liquid residuals, and many small waves. (The conditions under 
which these waves may be described by a small amplitude approximation are an 
open question.) Only in this regime does drain radius appear to be a sig- 
nificant parameter. The capillary dominated regime merges continuously into 
another which is characterized by few waves of large amplitude which repeat 
with small damping. This region is important, first, because it can be 
reached experimentally and, second, because low gravity phenomena occur on a 
larger scale than in the capillarity regime and more rapidly than in the 
draining dominated regime. 

For the fill levels considered here, high Weber number cases fall into this 
regime, which is characterized by a sheet of liquid left behind on the wall 
as liquid is withdrawn rapidly from the center of the tank. Fill level is 
an important parameter; with large fills there is more time for the contact 
point to fall and for liquid retained high on the wall to drain into the 
bulk of the liquid. The one example of this phenomenon encountered 
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in this study did not produce geysering but rather a falling, thinning wall 
sheet. The simulation produced a moving free surface that approaches the wall 
so closely that, in a corresponding experiment, the liquid would separate, 
leaving a sheet on the wall, and the bulk of the liquid would continue into 
the drain. This behavior violates a fundamental assumption in the model: 
the liquid remains connected. Geysering was provided for; but not separation. 
The only problem in including separation in the simulation program is an 
adequate physical description of the separation process. 

Because analysis was limited by the time needed to develop the program and 
learn how to use it to produce the simulations, all data files generated in 
the simulations have been stored on magnetic tape. Computer assisted re- 
analyses of existing data is possible. The frequent restart points can be 
used with the present or an augmented program to complete a simulation or to 
obtain additional analytic data for selected time intervals from a restart. 

During the course of this study a special descriptive terminology was adopted 
for referring to waves in the streamline end trajectories. A wave begins at 
the free surface when the trajectories near the contact point are at a maxi- 
mum displacement to the left; it reaches its maximum when the maximum dis- 
placement of the streamline ends to the right reaches the centerline; and 
it begins again when the maximum displacement to the left begins at the wall. 
Thus the beginning of a wave corresponds to a hang-up or slowdown at the wall, 
the maximum to a hang-up or slowdown at the centerline. 

The discussion of the data begins in the capillarity dominated regime and 
moves through the finite amplitude region at Bond number 0 to the development 
of wall sheets at high Weber number. This is followed by a detailed discus- 
sion of the thinning wall sheet case, with recommendations for continuation. 
The effect of Bond number is discussed next. 

The effect of fill level is considered last, because two of the simulations 
involved are those compared with experiments in the part which follows im- 
mediately. 

Rises at the Wall . The impulsive opening of the drain sets up instantaneously 
throughout the liquid, a velocity field directed toward the orifice. - Because 
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the velocity potential satisfies Laplace's equation, the initial streamlines 
and velocity field are perpendicular to the free surface (cf. Figure 7(c)). 
For the fill levels considered here, the initial streamlines depend princi- 
pally on Bond number, are curved near the free surface, are straight lines 
parallel to the axis in the body of the tank, and form a fan into the drain 
in the hemisphere (cf. fl] , pp. 6 3, 67 , and 95) fo r full initial streamline 
plots). Tracing a given initial streamline close to the wall backward from 
the body of the tank yields first a straight line segment parallel to the 
wall and then a curve that moves away from the wall to become perpendicular 
to the free surface. A large volume of liquid is trapped between the initial 
streamline, the free surface, and the wall. Because the streamlines move 
very slowly towards the wall and the free surface is hardly changing in this 
region, the volume Just described can be used as an instantaneous test volume 
(cf. Figures 7(c) through (h)). Because the lower part of the given stream- 
line is essentially fixed and parallel to the wall, liquid is flowing from 
this region into the drain at a constant rate. (Were a fixed test volume to 
be used, an insignificant outflow over the boundary within the tank would 
develop as the streamlines moved toward the wall. ) The radial component of 
the velocity at the free surface crams the upper part of the volume against 
the wall more rapidly than can be compensated for by the rate at which liquid 
is being withdrawn into the drain. Hence liquid near the contact point must 
climb the wall. This effect recurs repeatedly during wave motion whenever 
the streamlines near the contact point curve away from the wall (cf. Figures 
6(f), (g)> and (h)j and Figures 10(f) and (g)). Frequently, the net result 
is a marked slowing in the fall of the wall point. The initial rise at the 
wall, and the tendency to do so thereafter, depends principally upon the 
curvature of the free surface and the behavior of the streamlines near the 
contact point. Therefore, the occurrence does not depend upon the contact 
angle. However, the magnitude of the effect may be affected by the choice of 
contact angle. (A contact angle smaller than 5 °, possibly 0°, can be used in 
simulation. ) 

The Initial Transients . In response to the impulsive opening of the drain, 
the entire free surface falls, everywhere and in every case. At the contact 
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point for low Weber numbers, this response may be visible only in the initial 
velocity distribution; at higher Weber numbers the contact point may continue 
to fall for a time. The magnitude of the fall is so small that, even when 
it is visible in the time history of the wall point, it appears as an initial 
straight line segment. Both the magnitude and duration of the fall appear 
to increase with Weber number and to decrease with Bond number. The geometric 
tendency to rise at the contact point eventually overcomes the initial down- 
ward velocity. At low Weber numbers the rise is small and brief; at large 
Weber numbers it is larger and may continue to vapor ingestion. The maximum 
rise height and the time at which it occurs appear almost independent of fill 
level; both are decreasing functions of Bond number and increasing functions 
of Weber number. 

Meanwhile, the centerline point on the free surface continues its initial 
fall. Although the downward velocity is decreasing very slowly, the initial 
segment of the time history of the centerline point plots as a straight line. 
The "slope" of this initial segment appears to be characteristic of draining 
at the given Bond number. More precisely, all initial segments at a given 
Bond number are initially coincident in shape, independent of fill level; the 
higher and closer together the two Weber numbers in a comparison, the longer 
the coincidence lasts. A large neighborhood of the centerline point on the 
free surface appears to retain its initial shape as it falls. This neighbor- 
hood is defined as that region of a subsequent streamline plot in which both 
the free surface and the streamlines can be brou^it into coincidence with the 
original free surface and streamlines. As draining progresses, this neighbor- 
hood contracts toward the axis. The time at which this contraction reaches 
the axis appears to be of the order of the length of time that the initial 
segment of the centerline time history remains coincident with one for higher 
Weber number. The initial stretching of the free surface takes place between 
the wall and the neighborhood which still retains its initial shape. Modifi- 
cation of the program to output and process data quantities which are already 
computed in the course of the simulation would aid in understanding the 
physical processes behind these observations. 
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Effect of Drain Radius . Case 1 (0, .01, l/lO, 5 rr/3)> large drain radius, 
and Case 2 (0, .01, l/30, 5 tx / 3)> small drain radius, were chosen as a pair 
likely to show the effect of drain radius. There is a small effect; however, 
it occurs only near the drain. Hie streamline trajectories in Figures 4(a) 
and 5(a) are identical throughout most of the simulation and do not begin to 
diverge until after the centerline height, h , has fallen below . 3. Then the 

C 

trajectories in both cases begin to approach the axis, but those for the 
small drain do so more rapidly than those for the large drain. This is also 
evident in Figures 4(b) and 5(b). Although the behavior at the wall is 
identical, the centerline curve for the large drain near termination is 
flatter than that for the small drain. Figure 20 shows the patterns of 
streamlines at comparable times very low in the tank and higher in the tank. 

The plots in this figure have been carefully aligned so that a triangle may 
be used to verify the position of comparable points. Figure 20(a) shows that, 
very low in the tank, the streamline pattern for the larger drain lies to the 
right of that for the smaller drain. Figure 20(b) shows that the streamline 
pattern for the two drain sizes coincides except near the drain. The two 
panels on this page are characteristic of the two drain sizes in the following 
sense: in any simulation when the free surface is well above h c = .5, the 

lowest part of its streamline pattern will coincide with one of the two panels. 
This is a consequence of the draining nondimensionalization, which gives all 
drains the same outflow rate and allows them to differ only in the size of 
the orifice. 

Waves at Low Weber Number and Bottom Effects . Cases 1 and 2 show eight pure 
capillary waves, excited by draining, before the free surface begins to move 
toward the drain. Figures 5(c) and (d) show. Case 2 on an expanded scale with 
more detail of the trajectories near the wall. The amplitudes of the excur- 
sions of the streamline end trajectories toward the wall are alternately small 
and large. A similar pattern can be found in the centerline and wall height 
plots. 

Case 29 (5> .01, l/30, 5 tt/ 3) differs from Cases 1 and 2 in that the Bond 
number has been increased to 5* Figures 10(a) and (b) show 16 draining in- 
duced capillary- gravity waves with small amplitude occurring before vapor in- 

i, 
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gestion. Detail near the wall is provided by the expanded scale plots, 

Figures 10(c) and (d). In contrast with the alternating amplitudes of the 
pure capillary waves in Case 2, the amplitudes of the capillary-gravity 
waves in Case 29 steadily decrease and then increase as the base of the 
vibrating streamlines progresses into the hemisphere. 

When the free surface is high in the tank, the streamlines are relatively 
widely spread throughout the tank at the level at which they begin to vibrate. 
Figures 6(d), (e), and (f) are examples of this aspect of general streamline 
patterns. The feature that changes in the complete sequence of streamline 
plots through the period in which the increase in wave amplitude occurs is 
that their lower parts are progressively packed by the hemispherical bottom. 

This is shown in Figures 10(e) and (f). In Figure 10(e), near the maximum of 
the wave of largest amplitude low in the tank, the streamlines are confined 
by the wall and packed throughout the tank. At the beginning of the next 
wave. Figure 10(f), the streamlines remain packed in the center of the cross- 
section; but because the contact point is still on the vertical wall, the 
component of velocity along the free surface is choked off with the relaxation 
of the streamlines at the wall. 

When the contact point is on the hemispherical bottom, the streamlines do not 
relax much at the beginning of the wave, Figure 10(g), and consequently the 
strong velocity component along the free surface is maintained. This con- 
tinues until termination. Figure 10(h), and accounts in part for the low 
residual in this case. Although a considerable velocity component may de- 
velop along the free surface when h c is low In the tank, it may be cut off by 
wave action so long as the contact point remains on the wall as in Case 2, 

Figure 20 (a) . 

Cases 1, 2, and 29 are the only simulations in which a considerable portion 
of the calculations occurred in the region where the evaluation of the auxiliary 
potential is unsatisfactory. This fact may cause or contribute to numerical 
oscillations in the derivatives at the centerline, which are exhibited graph- 
ically by the muddy terminations of the two streamline trajectories adjacent 
to the axis in these three cases. Another factor may be the number of ver- 
tical mesh spacings used in the final stages of the calculation: k, 5, and 7 
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for Cases I, 2, and 29, respectively. Overall satisfaction with the endings 
of these three cases increases as the number of mesh spaces used increases. 

'Hie percentage change between original and final volume is much higher for 
Cases 1 and 2 (and presumably Case 29), . 276$, than the average of .02U$ for 
the remainder of the cases. Although these simulations take 6000 (Case 29) 
or nearly 8000 (Cases 1 and 2) time steps, the length of the run does not 
seem to be the cause of the loss of accuracy. 

Wave Dispersion, Nonlinear Effects . Case 5 (°, *1, l/30, 5 7r/3) is -an example 
of capillary wave action with all details clearly visible. Figures 6(a) and 
(b) have many as yet unexplained features. When the minor oscillations are 
smoothed out, it becomes clear that the second major wave reproduces the 
essential features of the first. The central part of the run including the 
second maximum was repeated with a smaller time step and the quivers in the 
trajectories were reproduced. This case was chosen as the example of the 
sequence of meshes used in a simulation (Figure 3) to show that there is no 
problem with the mesh. The streamline patterns in Figures 6(c) and (d) are 
typical of the shape at maximum wall rise and at the maximum of a wave. Com- 
parison of Figure 6(e) with Figure 7(c), which is an example of the common 
initial shape of the free surface for all cases with B=0 and volume 5 tt/Z, 
shows that the free surface can recover its initial shape in the course of a 
capillary wave. This occurs just after the maximum of the first wave at the 
end of a period of especially rapid fall at the centerline. Figure 6(f) 
shows the streamline adjacent to the wall bent over very far toward the center 
just after the wall point has started to rise toward its second maximum. It 
is just before and during this rise that unexpected behavior begins. Figure 
6(g) shows the same Btreamline again far to the left during the rise at the 
wall just before termination. Figure 6(h) shows that at termination, although 
there is a component of flow along the free surface toward the axis, the 
residual is large because the wave action is trapping much liquid high on the 
wall. 

Comparing Case 5 with Case 2 (Figure 5(c)) shows that the gross features of 
the former are replicated in the latter, but are somewhat smoothed. 
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Development of Wall Sheets . Case 8 (0, 1 , 1 / 30 , 5 77 / 3 ) shown in Figures 7 (a) 
and (b) is either a case in which the wave is very spread out, with none of 
the capillarity features of the previous Bond number zero cases, or a first 
example of a case with a wall sheet, a draining dominated case. Hie center- 
line initially falls at a rate characteristic of Bond number zero and then 
transits to a second region in which the fall is again approximately constant. 
Hie straight line segment in the wall history before the wall point begins to 
rise is actually an interval in which the wall point falls in response to the 
initial draining impulse before it begins to rise. Hie comparison of Figures 
7 (c) and (d) side by side may be difficult, overlaying shows that although 
the free surface in the latter has stretched at the wall, it still has its 
initial spherical shape at the center. Figures 7(e), (f), and (g) are a 
sequence showing the streamline shapes as the wall sheet develops and the 
centerline decelerates in response to a component of flow along the free sur- 
face. Figure 7(b) shows the streamline pattern just after the wall begins to 
fall. Note that at this time the center is well into the second interval of 
constant velocity. Figures 7(0 and (j) show the patterns at the beginning 
and end of the. final acceleration into the drain. 

Case 10 (0, 10, l/lO, 5 n/ 3 ), Figures 8 (a) and (b). Case 19 (0, 10, l/lO, 

8 77/3), Figures 9(a) and (b). Case 37 ( 5 , 10, l/lO, 5 tt fj>), Figures 13(a) 
and (b). Case 40 ( 5 , 50, l/lO, 5 77/3), Figures 14(a) and (b), and Case 49 
( 5 , 50, l/lO, 877/3), Figures 17(a) and (b), are all characterized by the de- 
velopment of a thick wall sheet, a long continued rise at the contact point, 
and an initial fall of the centerline at the rate appropriate to the Bond 
number. Hie two larger fill cases develop a second region of constant fall 
at the centerline; but the smaller fill cases go immediately from the initial 
fall into vapor ingestion. In Cases 37 and 49, at Bond number 5, the wall 
does fall, but very slowly. 

Hie final free surface on Figure 13(a) has a large salient in the wall sheet 
and has jumped into the drain. The final curve is the result of the predictor 
step following the time step shown in Figure 13(d)* The preceding time step 
is shown in Figure 13 (c). These figures show the top three lines of the mesh, 
a wave at the base of the wall sheet, and indicate a mesh problem just above 
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the drain. It is probable that there should be a small ware, as shown, and 
that it has been magnified to instability as the result of trying to pack too 
much mesh into the region above the drain. The liquid residual in Table I 
is computed from the free surface of Figure 13(d). Because the percentage 
change in volume is .0 6U$, the result is reliable. Moreover, examination of 
the centerline time history 6hovs that the rerunning of the end of this case 
would merely result in a miniscule change in the liquid residual and reduc- 
tion of the percentage change in volume to the general level. 

The Thinning Wall Sheet. Case h6 (5, 10, l/lO, 8r?/3), Figures 16(a) and (b), 
behaves initially as a draining dominated case in which the wall point is 
falling, eventually quite rapidly. The centerline point exhibits a new be- 
havior. It starts falling at the slope characteristic of Bond number 5, then 
transits to a second constant rate of fall, which also occurs for Case U9 , at^ 
finally transits to a third (new) constant rate for a shorter period Just 
before it starts to accelerate into the drain. The values for this case in 
Table I are for the terminal point of the doubtful segment shown in Figure 
16(b). When the terminal segment of the centerline time history of Case 37, 
Figure 13(b), is used to extrapolate the centerline history of Case U6, the 
resulting estimated volume is at most 3$ lower than the value given in Table 
I (see the following paragraph on Effect of Fill Level ). 

Comparing Figure 16(c), the streamlines and free surface shape when the com- 
putation is running smoothly, with Figure 16(d), the streamlines at the time 
of the last free surface plotted on Figure 16(a), shows that the wall has 
fallen, shortened, and thinned considerably. In the last set of streamlines 
the .02 streamline reaches up into the wall sheet. Thus the wall sheet is 
providing at least two percent of the flow out the drain and this is going 
through the neck in the wall sheet. The computer listing shows that the maxi- 
mum vertical velocity downward on the free surface is at the centerline and 
that there is a secondary maximum at the mesh point Just below the neck. The 
table below shows the widths of the neck, the local maximum vertical velocity 
(at the point on the free surface below the neck), the maximum width of the 
wall sheet, its length (above the neck), and Its volume as a percentage of 
the current volume in the tank. 
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Neck 

Wall Sheet 

Time 

Step 

Time 

Width 

Maximum 

Vertical 

Velocity 

Maximum 

Width 

Length 

Volume 

685 

1.9707 

.006 

-1.41 

.012 

.932 

2.55* 

695 

2.0004 

.004 

-1.66 

.012 

.941 

2.64* 

702 

2.0213 

.002 

-2.10 

.012 

.953 

2.66* 


The volume within the wall sheet decreases by about 3# during the period 
covered by the table; while the total volume decreases by about 7$. The 
minimum vertical velocity on the free surface outside the wall sheet, which 
is the easily available measure of the bulk velocity in the liquid, increases 
from -1.12 to -1.55 in this period. The vertical velocity distribution along 
the free surface above the neck is a uniform decrease in magnitude from the 
neck to the contact point where it changes from -1.20 to -1.23. Thus both 
the wall sheet and the bulk of the liquid are speeding up; but the latter 
more than the former. It appears that the inviscid model is predicting that 
liquid will separate and the wall sheet will be left flowing down the wall 
at a slower rate than the bulk of the liquid. 

This unexpected behavior was not included in the design of the program; how- 
ever, it is implicit in the mathematical model, for there is no reason why 
the dynamics of the problem should not force the free surface against the 
wall. Were the physics for breaking the domain into two pieces included in 
the model it would not be difficult to adapt the program to carry on with 
the lower piece. The fate of the wall sheet could also be determined. 

The thinning of the wall sheet brings into question the assumption that vis- 
cous effects are negligible. Ihe Justification of this assumption is based 
upon the fact that experimentally observed wall sheets have thicknesses much 
greater than those predicted by the analysis given in Appendix A of [l] . 
There an asymptotic expression for the thickness of a viscous wall sheet in 
a circular cylindrical tank is given as 

h Q = 1.35 (/^ 0 /a) 2/3 
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for zero gravity conditions. In this expression p is the viscosity, r is 

characteristic velocity for the central downward flow in the tank, and h 

o 

is the dimensionless thickness of the wall layer. Recasting this expression 
into the familiar dimensionless groups, Bond (B), Reynolds (R), and Weber (w) 
numbers, gives 

= 1.35 (V/R) 2 / 3 

as an equivalent expression. With the assumption of a Reynolds number of 
10,000, the values from Case 46, B=5 and Wo 10, give h Q the value .01. The 
values for the width of the neck given above indicate that the neglect of 
viscosity may no longer be strictly justified when the free surface approaches 
the wall. However, the neglect may be the practical way to solve the problem. 

Effect of Bond Number . The most dramatic change with Bond number is the 
transition at the larger fill level from the draining dominated Case 19 to 
Case 46 just discussed. At the lower fill level the corresponding Case 10 
and Case 37 have essentially the same properties. 

The second most dramatic change with Bond number is the Jump from Case 8, 
which is draining dominated, to Case 35(5 >1> l/30> 5 n /3)> Figures 12(a) and 
(b), which has a single wave of finite amplitude. The combination of drain- 
ing and gravitational forces suppresses practically all visible capillarity 
features in the oscillations of the streamline end trajectories in Figure 
12(a). Figure 12(b) shows an initial rate of fall at the centerline that is 
characteristic of all Bond number 5 cases. 

In the transition from Case 5 to Case 32 (5> . 1, l/30* 5 n /3)> Figures ll(a) 
and (b), the waves double In number, become relatively smooth, and definitely 
decrease in amplitude. When Case 32 is compared with Case 29 (Figure 10(c)), 
some common features emerge: an initial straight drop of the free surface 

which lasts longer near the axis than at the wall, a jump toward the axis of 
the streamline end trajectories near the axis just before the maximum of the 
first wave, and an overshoot of the trajectories near the axis on the re- 
covery from the first wave, lhese appear to be characteristic of capillary- 
gravity waves because they occur with very small magnitude in Case 35* 
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Effect of Fill Level. There are three comparisons between fill levels in 
the draining dominated regime. The striking feature of all three is that 
the wall height curve for the smaller fill completely overlays that for the 
larger. Thus behavior at the wall is independent of fill level so long as 
wave action does not occur. At the centerline, the initial segment of the 
smaller fill level overlays that of the larger but -breaks away at different 
times. When large scale free surface and streamline end trajectory plots 
for the smaller fill are overlaid on the plot for the larger a curious 
pattern of agreement is revealed. In each of the cases, the streamline end 
trajectories for the smaller fill are identical with those for the larger 
for a short distance below the free surface. The area of disagreement 
spreads toward the axis and wall but never reaches the latter. The height 
at which the area of disagreement reaches the axis for the lower fill level 
is about 1.3 for Cases 10 and 19, Figures 8(b) and 9(b), 1.5 for Cases 37 and 
46, Figures 13(b) and 16(b), and 1.6 for Cases 4o and 49, Figures 14(b) and 
17(b). For each Weber number, the lower end6 of the centerline time histor- 
ies in the figures just cited can be made to coincide. What is more, the 
lower ends of the centerline time histories for all cases with the same Weber 
number coincide and the height of coincidence increases with Weber number. 

This fact was used in estimating how close Case 46 was to vapor ingestion. 

Case 43 (5, 1, l/lO, 8 ^/3)> Figures 15(a) and (b), is a higher fill level 
analog of Case 35* Even with different scales, the strong resemblance between 
Figures 15(a) and 12(a) is apparent; not only does the latter resemble the 
initial segment of the former but it also appears to be the terminal segment. 
When Figure 12(b) is laid over the initial part of Figure 15(b), the curves 
for the height at the wall can be made to coincide and the initial part of the 
curves for the centerline height are identical, diverging only when the curve 
for Case 35 begins to accelerate toward the drain. When an attempt is made 
to match Figure 12(b) with the terminal segment of Figure 15(b), the match is 
remarkable. The following results for Case 43 emerge from the last compari- 
son: the wall is falling at the same rate in both of the nearly straight 

segments; the rate of fall of the centerline in the second straight segment 
is approximately the initial rate of fall; and the terminal behavior is only 
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sHghtly different. The principal difference between the two waves in Case 
43 is that both the wall and centerline curves are flatter in the transition 
from the second wave to terminal behavior than in the transition between 
waves. Thus it appears that the effect of increasing the fill and keeping 
the other parameters fixed is to replicate the wave behavior starting from 
the beginning. Consequently, the value of the liquid residuals may depend 
strongly upon the initial fill level. 


Comparison with Experimental Data 

Data from two experimental runs made at the NASA -Lewis Research Center were 
supplied for comparison with the simulations. The difference in initial 
conditions is the principal problem in making such comparisons. The simu- 
lations begin with an equilibrium interface with zero velocity. The Lewis 
Research Center experimental procedure calls for initiation of outflow when 
the interface centerline reaches its low point in its first pass through 
equilibrium [2] . A second problem is the difficulty in reading the point 
at which the interface meets the tank wall. And finally there is the unknown 
effect replacing the 0° contact angle by the 5° contact angle used in the 
simulations. 

Points from an experimental run with the parameter set (5.0, . 905, .10, 5 77/3) 
are compared in Figure 18 with the computed centerline and wall height curves 
for Case 35 (5> 1> l/30, 5 n / 3) • The experimental points lie very close to 
the computed centerline heights; there is agreement in level and trend between 
the experimental and computed wall heights. At this Bond number, the differ- 
ence in' drain radius should have no effect; however, the experiment with the 
smaller Weber number should have a longer drain time; and this is observed. 

At the centerline the initial experimental point is low with respect to the 
simulation and it initially falls at a rate less than the initial rate charac- 
teristic of Bond number 5» This suggests a small upward velocity at the 
centerline. At the wall the experimental point is low and initially falls. 
This suggests a downward initial velocity at the wall. The divergence in the 
initial conditions accounts for the flattening of both experimental curves 
with respect to the computed curves. 
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Although it has not been specifically brou^it out, experience in following 
waves in streamline end trajectories across the tank suggests that the ex- 
cess experimental velocity at the wall travels across the tank to the center- 
line and back to the wall again. When all the excess velocity finally leaves 
the wall, it falls at the rate characteristic of the basic problem. The ex- 
cess velocity reaches the centerline just before the maximum of the wave 
occurs and accounts for the increase in the experimental rate of fall. The 
excess velocity returning to the wall again contributes to the divergence of 
the rates of fall at the wall just before vapor ingestion. Examination of 
Figure 18 suggests that the half period of the wave of excess velocity is 
about .5; from Table H, the value of the period for Case U 3 is . 9 6 . 

In spite of the difficulty in reading wall points, there is a coincidence of 
the experimental and computed curves during the time that the wall is falling 
at the rate characteristic of both Case 35 and Case 4 3 . After the effect of 
the small initial velocity at the centerline dies out, the experimental center- 
line point does fall at approximately the initial rate characteristic of Bond 
number 5 and follows the computed curve through the transition into the maxi- 
mum of the wave very well. Although the experimental points are lower, both 
curves fall at similar rates. Ihe chief discrepancy between experiment and 
simulation is that the latter predicts a slowing down of the wall point just 
before vapor ingestion, whereas the experimental point continues to fall as 
a result of the wave of excess velocity. The larger contact angle in the 
simulation may also contribute to the difference. The agreement in this com- 
parison is excellent. 

Points from an experimental run with the parameter set (5.0, .9^5, .10, 8 ff/ 3 ) 
are compared in Figure 19 with the computed centerline and wall height curves 
for Case ^3 ( 5 , 1, l/lO, 8 tf/ 3 ). Again the experimental points at the center- 
line lie close to the computed curve. However, the experimental determina- 
tion of the wall heights only resembles the computed curve in general slope 
and lies so much below it that to use the experimental wall heights in con- 
structing a conjectured initial velocity to explain the divergences from the 
good agreement at the centerline would be unrealistic. The initial centerline 
point is above the computed point. However, both fall at the rate character- 
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istic of Bond number 5 for approximately the same length of time. Thus it 
would appear that the initial experimental interface had the shape of the Bond 
number 5 equilibrium free surface and a very small residual velocity over a 
very large part of the tank. Ihe first experimental wave begins about the 
same time as the computed one reaches its maximum perhaps a bit earlier, and 
is of much snaller amplitude. Both sets of data return to an approximation 
of the characteristic rate of fall. Because it is expected, an experimental 
second wave of small amplitude can be detected at the appropriate time. The 
smaller experimental Weber number again produces a slightly longer drain time. 
Although the agreement is not so striking as in the previous comparison, it 
is good, particularly in the initial centerline detail. 

Liquid Residuals 

The quantity of liquid left in a tank after vapor ingestion occurs has con- 
siderable impact on the design of orbital transfer systems. Residual weight 
must be accepted as a weight penalty or traded with the weight and other dis- 
advantages of alternate means that reduce residuals below values for a clean 
tank. Obtaining data on residuals for low gravity draining of hemispherically 
bottomed, cylindrical tanks was a major objective of this investigation. The 
data in Table I and the plots in Figures 21 - 23 are useful for design 
studies and meet this objective, Hie following paragraphs are a discussion 
of the relative significance of drain size, initial fill volume, Bond number, 
and Weber number on the residual volume. Generally, the numerically deter- 
mined values of residual volume are somewhat larger than those obtained ex- 
perimentally [2] ; the differences are attributed, primarily to noneauilibriuir 
free surfaces in the experiments at the start of draining and, secondarily to a 
larger contact angle (5°) used in the calculated cases. The numerically de- 
termined values are conservative, however. 

Drain Size . Observations regarding the general influence of drain size on the 
draining solutions are discussed earlier in this section; here specific re- 
sults are presented regarding the effect of drain size on the residual volume. 

From the calculations of [l] , it was concluded that Weber number must be 
much less than one or that for a given Weber number the Bond number must be 
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much greater than one to cause appreciable parametric dependence of- residual 
volume on drain size. Cases 1 and 2, in which W = .01 and r Q = l/lO and 1/30, 
respectively, were selected to investigate this parameter. (The same values 
hut with B=5 would have shown a greater effect, but data for B=0 were con- 
sidered particularly desirable.) Review of Figures 4(a), 4(b), 5(a), and 
5(b) do not show appreciable differences in draining features between Cases 
1 and 2 up to vapor ingestion; this is particularly evident in Figures 20(a) 
and 20(b) in which streamline patterns are compared. Had a better numerical 
ending been possible for Case 1 (cf. Figures 4(a) and 4(b)), a larger numeri- 
cal difference between values of V would have been found. The values of V 

r r 

for Cases 1 and 2, which are given in Table I, differ by 1.5$* Thus, one 
must reduce below the values for Case 1 and 2 to obtain appreciable varia- 
tion of residuals with drain size. 

Initial Fill Volume. The initial fill level, h , defined as the initial 

— — c 

height of the equilibrium free surface at the tank axis, was chosen as a 

parameter in the survey of draining solutions reported in [ 1 ] . The initial 

fill volume, V , is used herein instead. There are advantages in selecting 

either one, but the data for given in [ 1 ] suggest that a survey at fixed 

values of V offers a more consistent and useful grouping of residual volume 
o 

data when considering the major parameters of Bond number and Weber number. 
Also, in applications one is interested in a draining efficiency, 1 - V ' /V , 
so that the best values of Bond and Weber number can be combined to minimize 
this quantity. Two values of V q were selected, namely, 5 77/3 and 8 nf 3 (cf. 
Table i). For fixed values of B, W, and t q , Figure 8 (Case 10) is compared 
with Figure 9 (Case 19 ); Figure 13 (Case 37) with Figure 16 (Case 46); and 
Figure l4 (Case 4o) with Figure 17 (Case 49). Except for Cases 37 and 46, 
the residual volume increases as V q increases; evidently that the longer wall 
layers for the deeper cases cause the larger residual volume. Generally, as 
W and V increase residuals increase, but one must use caution because for 
some intermediate Weber numbers and nonzero Bond number, such as in Cases 37 
and 46, the residual volume decreased instead. For these two cases the wall 
point starts to move downward quite rapidly after the characteristic initial 
delay, and the resulting wave action on the free surface extends the time for 
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vapor ingestion and thus gives lower residuals. The conclusion that residuals 
increase with increasing initial fill volume is, therefore, conditioned by 
saying that the Weber number must be so- large that the wall point never has 
enough time to get started moving downward at an appreciable rate before vapor 
ingestion occurs. 

The most remarkable results obtained from varying initial fill volume are for 
Cases 35 and 43 (cf. Figures 12(a) and 15(a)). Case 35, where = 5 n / 3, 
appears to duplicate Case 43 once the free surface in the latter reaches a 
value of z = 1.7 to 1.8 at the centerline. From Table I the values of residual 
volume differ very little suggesting that for W=1 and B=5, the transition has 
been made (cf. Figure 2l). This agreement may be coincidental for the two 
cases because the residuals will be altered if one has vapor ingestion in 
different parts of a wave. As the Weber number decreases, this effect should 
diminish and the dependence on fill level should disappear. Unfortunately, 
deep cases for B=0 are not available to define the transition, but comparing 
Case 5, in which W = 0. 1, and Case 8, in which W = 1, suggests that the 
transition likely occurs at or near W = 0,1 for E=»0 because the appearance 
of the wave action in Figure 6(a) is so similar to that in Figure 12(a) for 
B = 5. For Weber numbers much less than one, ample time is available to allow 
the free surface to oscillate freely, producing the patterns seen in Cases 
1, 2, and 29, for example, and eliminate any dependence on V . 

Weber Number, Bond Number, and W/(l+B ). Several cases for B = 0 and 5 were 
computed to define the dependence of residual volume on Weber number. The 
data reported in Table I and plotted in Figure 21 span three and one-half 
decades of Weber number. In retrospect, additional computations for even 
smaller values of W and a complete survey of W at a larger Bond number would 
have been desirable. However, the existing data, coupled with available ex- 
perimental measurements for very large Bond number [2] , suffice to define 
the pattern of the residual volume curves (Figures 21 and 22). 

As a function of W, the computed data for residual volume divide into three 
regimes: one, in which draining forces dominate (W large); two, in which 

capillary forces dominate (W small); and three, the transition between high 
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and low W (V of order one). Boundaries between the three regimes are not 

sharply drawn; but the asymptotic behavior in the large W and small W zones 

is evident from the S-shaped curves that result (see plots of 7 vs. W in 

r 

Figure 2l) and in the remarkably different character of the simulation plots 
(Figures 5 (a) - 17 (a)). The appearance of these regimes stems from the time 
scaling of a given problem. In one extreme, draining is so rapid that liq- 
uid is withdrawn only from the central part of the tank leaving wall layers 
and a large residual; in the other extreme, capillary and gravity forces 
have ample time to restore the free surface by forcing draining from the wall 
region as well as from the central part of the tank. The latter idea is amply 
demonstrated by the induced wave motion and oscillatory streamlines in 
Figure 10 (e) - (h). 

For engineering calculations in which an estimate of residual volume is re- 
quired, the graphs in Figure 22 are recommended. Here residual volume is re- 
lated to the parameter, W/(l+B), suggested in [l] as the simplest combination 
of W and B that produces the proper limits for B tending to zero and infinity. 

This parameter dramatically reduces the parametric dependence of V on B and 

r 

permits comparison with high Bond number data. Some parametric spread of the 
plots in Figure 22 with Bond number is expected because of variations due to 
low gravity free surface shapes which cannot be normalized with this simple 
relationship. 

All the cases computed in this investigation, as well as in [l],are shown on 
the plot of t vs. w/(l+B) in Figure 23. Using the formula V =3((v /n) -t )/2 
one can back out estimates for residuals regardless of initial fill volume^ 

V q ,or fill level, h c . Earlier work for high gravity draining resulted in 
simple expressions for computing a critical height defined as the liquid 
level at the incipience of vapor ingestion (see [2], [17] , and [l8 ] ). This con- 
cept for large B is very useful for predicting residuals because the shape of 
the free surface is essentially flat; for low-to-intermediate values of B 
thiB approach is not accurate and the substitution of estimates based on 
Figures 22 and 23 is recommended. 
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Wave Motion 


Cases with W/(l+B) < 1 (Cases 1, 2, 5, 29, 32, 35, and 1+3) developed axi- 
symmetric slosh at the free surface. For these cases the draining time is 
long enough to observe one or more periods of the wave motion. From the list- 
ings of the solutions, tabulations of the periods were made for those cycles 
not immediately altered by the rapid acceleration of the free surface into 

the drain. The values were averaged and the results tabulated as A , in 

di. 

Table II. Further, these wave periods were plotted against (W/(l+B)) 2 , the 
ratio of the sloshing and draining time scales (Figure 24). The straight 
lines in Figure 24 and their 45° slope suggest that the periods measured in 
the time scale of the draining nondimensionalization should be rescaled in 

n JL 

terms of the sloshing time unit, ( p&^/a (l+B)) 2 , introduced in [ 17 ] and 

discussed in Numerical Analysis as an alternate to the draining time unit, 

3 / * 

rt a-yQ. This was done and the results tabulated as A in Table II. The 

. * 

values of A are constant for the two Bond numbers; dispersions are attributed 

to the crude calculations used to obtain A . . 

d 

These results suggest that the periods of draining induced oscillations 

cluster about a constant value for each Bond number if Weber number is less 

than or equal to one. At higher values of Weber number the waves are still 

present but not easily observed because of short drain times compared to the 

period A^. For large Weber numbers and very deep cases which would prolong 

the drain time, one could expect large amplitude slosh with potential geyser- 

* 

ing at the centerline. In this case the dynamics would be nonlinear and A 
would likely exhibit considerable deviation from the values in Table II. 

In [ 1; Figure 12, pp. 95-99] , wave motion is simulated for B = 19*32, 

W = 3.41, and h =2.0. The value of A A is estimated to fall in the range 

1.2 - 1.3 and, as a result, A in 2.9 - 3 * 2 . Used in conjunction with the 

* 

values for B = 0 and 5 in Table II, A generally increases with Bond number. 

The investigation of Saad and Oliver (lo] demonstrated the existence of 
draining-initiated waves by a linear analysis in which -016 low-gravity 
liquid- vapor interface was assumed flat ( 90 ° - contact angle or a sufficiently 
large Bond number). Their work did not require the velocity-squared term in 
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the Bernoulli equation to obtain linearized wave action. For w/(l+B) small 
enough and small contact angle, an analysis similar to that of Saad and Oliver 
would provide results similar to those given herein because the simulated 
wave amplitudes are small. For W/(l+B) of order one or greater and with 
arbitrarily deep initial fills, a full simulation iB required to follow the 
finite amplitude motionj but additional simulations with the present computer 
program would be needed to define completely the transition regime between 
these two extremes. 
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SUMMARY OF RESULTS 


The basic results of the Numerical Simulation of Low Gravity Draining study 
are fourteen simulations and a computer program using new techniques. 

(1) Fourteen nondimensional simulations cover three and a half decades 
in Weber number, W a .01 through 50, two Bond numbers, B = 0 and 5, 
two drain radii, l/30 and l/lO of the tank radius, and two initial 
volumes, 5^/3 and 8 tt / 3 . 

(2) The computer program solves the difficult, but well posed, initial 
boundary-value problem over a wide range of parameters and forms of 
solutions. The continuity of the solutions as a function of the 
parameters is a fundamental tool used in the analysis to bridge gaps 
between widely spaced data points. 

( 3 ) Two simulations are compared with two similar NASA-Lewis Research 
Center experiments. Agreement is excellent. 

(4) The range of simulation in Weber number extends that of experiment 
by two decades. 

( 5 ) A principal goal of the study was to determine the liquid residual, 
V^, for a wide parameter range. Together with values from earlier 
Lockheed computations, is plotted against W and W/(l+B). The 
latter is the simplest parameter that connects somewhat similar 
cases and occurs in the transformation between the time scales used 
for nondimensionalizing draining and sloshing problems. 

( 6 ) For design studies, the recommended way to determine is to use 
the equation = 3 (( v 0 /ir) - t^)/ 2 , where is initial non- 
dimensional volume and t vi is vapor ingestion time given by the 
plot of against W/(l+B). 

( 7 ) The S-shaped curves of V against W suggest three regimes: (i) a 

r 

capillary regime, characterized by low Weber numbers, low residuals, 
and many small waves, in which capillary and gravitational forces 
have time to drain liquid repeatedly away from the wall; (ii) a 
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draining regime, characterized by high Weber number, large residuals, 
and the development of a wall sheet, in which the draining simply 
draws liquid from the center of the tank with the wall sheet left 
behind; and (iii) a transition region between the two extremes, 
where there are few waves of large amplitude. 

(8) Drain radius has a significant influence on residual volume for 
Weber numbers much less than one or, at a fixed Weber number, for 
Bond numbers much greater than one. 

( 9 ) The periods of the waves become a constant, characteristic of the 
Bond number, when normalized to the sloshing time scale. 

(10) Fill level is a very important parameter in the draining dominated 
regime. What happens at very high fill levels when there is time 
for the wall sheet to reorient is an open question. 

(11) Increasing the Bond number may (i) introduce more waves of smaller 
amplitude in the capillary regime, (ii) change the nature of the 
problem from one in the draining regime to one with waves, or 
(iii) simply reduce the time scale in the draining regime, that is, 
permit the wall sheet to thin more. 

(12) Control of mesh spacing along the stretching and contracting free 
surface was achieved by the method of variable trajectories in which 
a small auxiliary kinematic velocity is added to the Lagrangian 
velocity components at each time step to insure that the computa- 
tional points maintain a desired distribution. The Bernoulli equa- 
tion at the computational points is modified to account for the 
change in direction. This new method was developed specifically 
for this problem. 

( 13 ) The success of these simulations has demonstrated the utility of 
the method of moving triangular meshes, in which parts of the mesh 
are revised periodically and parts of mesh rows and columns are 
dropped and added to accommodate the motion of the free surface. 
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(l4) A simulation for B = 5 , W = 10, and an initial fill of 8 iff 3 pre- 
dicts on the basis of the inviscid, irrotational model that the liquid 
in a very thin falling wall sheet may separate and be left on the wall, 
while the bulk of the liquid continues down the drain. This simulation 
shows the power of the program; it continues the solution until the 
neck is .002 of a tank radius in width. It also suggests that there is 
a region near the boundary of the transition regime where the neglect 
of viscosity may be questioned. 
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Figure 1. Definition of geometry and coordinate system for the draining problem. 
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(a) Physical mesh at t=.8267. 


Figure 2. Triangular mesh generated for Case 37; 
(B,W,r o ,V Q )=(5, 10 , 1/10,5 tt/3). 
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(a) Initial mesh used for B=0 cases. 


Figure 3* Evolution of computing mesh for Case 5; 
(B,W,r o ,V Q ) = (0,. 1,1/10,5 tt/3). 
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(b) Mesh at t=.4^39. 


Figure 3 (cont). Evolution of computing mesh for Case 5; 


(B,W,r o ,V o )=(0,. 1,1/10,5 tt/3). 
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(c) Mesh at t=.886o. 


Figure 3 (cont). Evolution of computing mesh for Case 5 
( B > w ,r o ,V o )=( 0 ,.l, l/lO, 5 n/z). 
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Figure 4. Draining characteristics d 
(B,W,r o ,V o )=(0,.01,l/l0,5 
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centerline and wall. 

Figure ^ (cont). Draining characteristics for Case 1; 

(B,W,r o ,V o )=(0,.01, 1/10,5 tt/3). 
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(a) Free surfaces and trajectories of free 
surface and streamline intersections. 


Figure 5. Draining characteristics for Case 2; 
(B,W,r . V )=(0, .01, 1/30, 5 7T /3). 
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centerline and wall. 

Figure 5 (cont). Draining characteristics for Case 2; 

(B,W,r o , v o )=( 0 , .01, 1/30 , 5 7r/3). 
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(c) Expanded scale. 


Figure 5 (cont). Draining characteristics for Case 

(B,W,r o ,V o )=(0,.01,l/30,57r/3). 
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(d) Expanded scale. 


Figure 5 (cont). Draining characteristics for Case 2 ; 

(B, W, r 0 , V q )=(0, . 01, 1/30, 5 7T /3) . 
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(a) Free surfaces and trajectories of free 
surface and streamline intersections. 

Figure 6. Draining characteristics for Case 5; 
(B,W,r o ,V Q ) = (0,. 1,1/30, 5 tt/3). 
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Figure 6 (cont). 


Draining characteristics for Case 5; 
(B,W,r ,V )=(0,. 1,1/30, /3). 
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(c) Streamlines at t=.1201 

showing maximum rise at wall. 


(d) Streamlines at t=. 1+275 near maximum 
of first major wave showing large 
excursion of the .98 streamline and 
continued acceleration along the 
remainder of free surface. 


Figure 6 (cont). Draining characteristics for Case 5; 

(B,W,r o ,V o )=(0,. 1,1/30, 5 tt/3). 
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.0 .5 1.0 

(e) Streamlines at t=,5393 

resembling initial conditions 
(cf. Figure 7(c)) except at 
the wall. 


.0 .5 1.0 

(f) Streamlines at t=.6l38 near 
beginning of second wave 
showing the .02 streamline with 
curvature exceeding its initial 
curvature (cf. Figure 7(c)). 


Figure 6 (cont). 


Draining characteristics for Case 5; 
(B,W,r o ,V o )=(0,.l, 1/30,5 jt/3). 
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(g) Streamlines at t=l.l 831 near possible 
beginning of third major wave showing 
very curved streamlines adjacent to wall 
(cf. Figure 6(f)). 

Figure 6 (cont). Draining characteristics for Case 5; 

(B, W,r o ,V o )=( 0 ,. 1,1/30, 5 »r/ 3 ). 
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Figure 6 (cont). Draining characteristics for Case 5; 

(B,W,r o ,V o )=(0,.l, 1/30,5 tt/3). 
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Figure 7» Draining characteristics for Case 8; 
(B,w,r Q ,V o ) = (0,1,1/30,5 h-/3). 
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(c) Initial streamlines and 
spherical free surface. 


^ ' \ "■ 
\ * 


\ ' \ • 
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(d) Streamlines at t=.1086, when 
free surface is dropping at 
initial rate (cf. Figure 7(b)) 
showing that center is still 
spherical. 


Figure 7 (cont). Draining characteristics for Case 8; 

( B ,V,r o ,V o )=(°,i, 1/30,5 tt/3). 
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(e) Streamlines at t=.252k 
as center begins to 
decelerate. 


(f) Streamlines at t=. 3962 
showing developing 
velocity component along 
free surface. 


Figure J (cont). 


Draining characteristics for C-se 8; 
( B ’ W > r 0 ,V o )=(0, 1,1/30, 5 tt/3). 
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(g) Streamlines at t=.5l8U 
near maximum of wave. 
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(h) Streamlines at t=. 633 ^ just 
after wall begins to fall. 

Center line velocity is again 
nearly constant (cf. Figure 7(b)), 


Figure 7 (cont). Draining characteristics for Case 8; 

(S, W, r Q , V q )= (0, 1, 1/30, 5 7T /3) - 
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(i) Streamlines at t=.83^8 as (j) Terminal streamline pattern 

center line begins to at t=.9929. 

accelerate (cf. Figure 7(b)). 

Figure 7 (cont). Draining characteristics for Case 8; 

(B»W,r o ,V o )=(0, 1,1/30, 5 tt/3). 
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(a) Free surfaces and trajectories of free 
surface and streamline intersections. 
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Figure 


8. Draining characteristics for Case 10; 
(B,W,r o ,V o )=(0,10, 1/10,5 tt/3). 
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(b) Histories of free surface points 
at centerline and wall. 


Figure 8 (cont). 


Draining characteristics for Case 10; 
(B,W,r o ,V Q )=(0,10, 1/10,5 tt/3). 
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(a) Free surfaces and trajectories of free 
surface and streamline intersections. 

Figure 9. Draining characteristics for Case 19 ; 

(B, W, r 0 , V q )= (0, 10, 1/10, 8 7 T / 3 ) . 
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(b) Histories of free surface points at 
centerline and wall. 


Figure 9 (cont). Draining characteristics for Case 19; 

(B,V, r o ,V o )=(°, 10 , 1/10,8 tt/3). 
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(b) Histories of free surface points 
at centerline and wall. 

Figure 10 (cont). Draining characteristics for Case 29; 

(B,W,r o , V q )=(5,. 01, 1/30,5 tt/3). 
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(c) Expanded scale. 


Figure 10 (cont). 


Draining characteristics for Case 29 

(B,W,r o ,V o )=(5,. oi, 1/30, 5 tt/3). 


LOCKHEED MISSILES & SPACE COMPANY. INC. 



STEP TIME 


4640. 1 . 2509 
4800. 1.2946 

4960. I . 3383 
5100. 1 . 3765 

5260. 1.4201 
5400. 1 . 4583 

5560. 1.4935 
5720. 1.5271 

5850. 1.5544 
5990. 1 . 5775 

6140. 1.5912 
6153. 1.5920 




(d) Expanded scale. 
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Figure 10 (cont). Draining characteristics for Case 29 

(B, W, r o , V q )= ( 5, . 01, 1/30, 5 tt/3) . 
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(e) Streamlines at t=l. 1718 near maximum of 
wave with amplitude enhanced by approach 
to hemispherical bottom. 


Figure 10 (cont). Draining characteristics for Case 29; 

(B,w,r o ,V o )=(5,.01,l/30,5 ?r/3). 
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(f) Streamlines at t=1.22 6 k near beginning 
of wave with amplitude enhanced by ap- 
proach to hemispherical bottom. 


Figure 10 (cont). Draining characteristics for Case 29 

(B,W,r o ,V o )=(5,.01,l/30,5 n/3). 
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(g) Streamlines at t=l. U 17 U showing beginning 
of wave with diminished amplitude in 
hemisphere. 


Figure 10 (cont). Draining characteristics for Case 29 

(B,W,r o ,V o )=(5,.01, 1/30,5 7r/3). 
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(a) Free surfaces and trajectories of free 
surfaces and streamline intersections. 


figure 11. Draining characteristics for Case 32; 
(B, W, r Q , V o ) = ( 5 , . 1, 1/ 30, 5 ?r/3) . 
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(a) Free surfaces and trajectories of free 
surfaces and streamline intersections. 



Figure 12. Draining characteristics for Case 35; 
(B, V, r 0 , V> )«(5, 1 , 1/30, 5 rr/3) . 
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(b) Histories of free surface points at 
centerline and vail. 

Figure 12 (cont). Draining characteristics for Case 35 ; 

(B, v,t o ,V o M 5, 1,1/30, 5 rr/3). 
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5 

(a) Free surfaces and trajectories of free 
surface and streamline intersections. 


Figure 13 . Draining characteristics for Case 37; 
(B,W,r o ,V o )=(5,10,l/l0,5 tt/3). 


LOCKHEED MISSILES 8c SPACE COMPANY. INC. 





(b) Histories of free surface points at 
centerline and wall. 


Figure 13 (cont). 


Draining characteristics for Case 37; 

( B >w,r o ,V o M5, 10 , 1/10,5 
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(c) Mesh and free surface at t=. 8712 
showing wave at base of wall sheet. 

Figure 13 (cont). Draining characteristics for Case 37; 

(B,W,r o ,V o )=(5, 10, 1/10,5 n/3). 
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(d) Mesh and free surface at t=. 87^-2 
just before vapor ingestion. 

Figure 13 (cont). Draining characteristics for Case 37; 

(B,W,r o , V q )=(5, 10, 1/10, 5 tt/3). 
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(a) Free surfaces and trajectories of free 
surfaces and streamline intersections. 

Figure 14. Draining characteristics for Case hO; 

(B,W,r ,V )=(5, 50, 1/10,5 tt/3). 

106 o o 


LOCKHEED MISSILES & SPACE COMPANY. INC 






(b) Histories of free surface points at 
centerline and vail. 

Figure 14 (cont). Draining characteristics for Case 40; 

( B ,W,r o ,V o M5,5 0 ,l/l°,5 n-/3). 
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(a) Free surfaces and trajectories of free 
surface and streamline intersections. 

Figure 15. Draining characteristics for Case ^3; 
(B, W, r Q , V q )=( 5, 1> 1/10, 8 tt/3) . 
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(b) Histories of free surface points at 
centerline and wall. 

Figure 15 (cont). Draining characteristics for Case L3; 

(B,W,r o ,V o )=(5, 1,1/10 , 8 tt /3). 
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(a) Free surfaces and trajectories of free 
surface and streamline intersections. 

Figure 1 6 . Draining characteristics for Case 46; 
(B,W,r ,V )=(5,10, 1/10,8 tt/3). 
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Figure 16 (cont). Draining characteristics for Case U6; 

(B,W,r o ,V o )=(5,i0, 1/10,8 tt/ 3 ). 
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(d) Streamlines at t=1.9hll shoving 
necking at bottom of vail sheet. 


Figure 1 6 (cont). Draining characteristics for Case k6 

(B, W, r Q , V o ) = ( 5 , 10, 1/10 , 8 n /3) . 
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(a) Free surfaces and trajectories of free 
surface and streamline intersections. 

Figure IT. Draining characteristics for Case k9; 
(B,W,r o ,V o )=(5, 50, 1/10,8 ^/3). 
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(b) Histories of free surface points 
at centerline and wall. 

Figure 17 (cont). Draining characteristics for Case 49; 

(B, w , r 0 , V Q )= ( 5, 50, 1/10 , 8 tt/3) . 
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Figure 18. Comparison of computed and experimental data for Case 35; 

(B ' W > r o'V = (5,l,l/30,5rr/3). 

NASA parameters are (5.0, .905, .10, 1.67/0. 
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C43 B = 5.0 W = 1.00 VOL = 2.67 DR = .1000 

* NASA B = 5.0 W = .945 VOL = 2.67 DR = .1000 


WALL 
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Figure 19- Comparison of computed and exDerimental data for Case ^ 3 ; 
(B,W,r 0 ,V o ) = (5,1,1/10,8 tt/3). 

NASA parameters are ( 5.0, .9^5, .10, 2 . 6777 ). 
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(a) Pattern for r =l/lO lies to right 

of that for r°=l/30 when free surface 
is near drain? 


Figure 20. Comparison of streamline patterns for large 
and small drain radii. 
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STEP TIME 
6000. 1.2239 



(b) Patterns coincide except near drain 
when free surface is well above drain. 

Figure 20 (cont). Comparison of streamline patterns for large 

and small drain radii. 
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TABLE I SUMMARY OF COMPUTED DATA 

(All quantities are dimensionless) 


Case 

number 

Problem parameters 
(B, W, V V o ) a 

Residual 

volume (V ) 
r 

Vapor ingestion 
time (t^) 

Wall height at 
vapor ingestion 

1 

« 

(o, .oi, i/io, 5*73) 

.28752 

1. 4789 

1.131 

2 

(0,. 01, 1/30, 5 rr/3) 

.29203 

1.4761 

1.136 

5 

(0,.l, 1/30,5 ?r/3) 

. . 571^8 

1.2888 

I.433 

8 

(0,1,1/30,5^3) 

1.00629 

• 9958 

2.653 

10 

(o,io, 1/10,5*73) 

- 1. 45462 

.6967 

2.676 

19 

(o,io, l/lO,8*/3) 

1.67175 

1.5523 

3.856 

29 

(5,. 01, 1/30,5*73) 

.11374 

1.5919 

.614 

32 

(5,.l, 1/30, 5 */3) 

.24067 

1.5097 

.798 

35 

(5,1, 1/30, 5 * 73 ) 

.70750 

1. 1989 

1. 356 

37 

(5,10, 1/10,5*73) 

1.18648 

.8742 

2.308 

40 

(5, 50, 1/10,5*73) 

1.25088 

.8326 

2.430 

43 

( 5,1,1/10,8 */3) 

.70353 

2.2053 

1. 346 

46 

(5,10,1/10,8 jr/3) 

.96555 

2.0213 

2.360 

49 

(5, 50-, l/l0,8 */3) 

1.44102 

1.7063 

3. 399 


a (B, W, r , V ) = (Bond number, Weber number, drain radius, initial liquid 
0 ° volume) 
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WEBER NUMBER (w) 

Figure 21. Computed values of residual volume for low gravity draining of 
heraispherically bottomed, cylindrical tanks. (Data for initial 
fillings of 5 7 t/ 3> dashed symbols for initial fillings of 8 7l/3. ) 











Figure 22. Computed values of residual volumes for low gravity draining 
of hemispheric ally bottomed, cylindrical tanks. (Data for 






DRAINING PARAMETER (W/(l+B)) 

Figure 23. Computed vapor ingestion times for low gravity draining of hemispherically 
"bottomed, cylindrical tanks. (Data for h = 1.0, 2.0, 2.6652, 2.7l85> and 
3.0 taken from [l ] ^ Table I.) 










0*1 



RATIO OF DRAINING AND WAVE MOTION TIME SCALES (w/(l+B )) 3 

Figure 2k. Period of free surface waves produced by draining at low 
Weber numbers. 







TABLE II ESTIMATED WAVL PERIODS 

(All quantities are dimensionless) 


Case 

number 

Weber 
number (W) 

Bond 

number (B) 

(W/(l+B))* 3 

Wave 

period (A,) 
d 

Wave „ 

period ( A ) 

1 

.01 

0 

.1 

.18 

1.8 

2 

.01 

0 

.1 

.18 

1.8 

5 

.1 

0 

.316 

.59 

1.9 

29 

.01 

5 

.0408 

.094 ’ 

2.3 

32 

.10 

5 

.1291 

.30 

2.4 

k 3 

1. 

5 

. 

fr 

O 

03 

.96 

2.4 


Ratio of sloshing and draining time scales. 

b A* = A d ((l+B)/W)2. 
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APPENDIX A 


SYMBOL LIST 


Quantities are nondimensional unless otherwise designated. When appropriate, 
the relation between a dimensionless variable and the physical dimensional one, 
which is topped by a bar, is given. Underlined variables are generally vectors. 
Spatial and temporal variables appearing as subscripts to R, Z, N, S, and V 
indicate partial differentiation. 


English Alphabet 
a 


B 

b 

c 

C(t) 

DR 

f 

f 

w 

F 

g 

h 


-Dimensional radius of container. The characteristic 
length with respect to which the variables are made 
nondimensional. 

-Arbitrary constant. 

-Bond number = pga /o. 

-Arbitrary constant. 

-Additive constant. 

-Additive function of time. 

-r on plotted output. 

-Height of point on instantaneous free surface meridian = f/a . 

-Value of f on container wall = f /a . 

w ; 

-Height of point on equilibrium free surface meridian = F/a . 
-Dimensional acceleration due to gravity. 

-Height of equilibrium free surface at the container 

axis = h /a . 
c' 

-Height of instantaneous free surface at container axis = h /a . 
c d 
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w 


W, VI 


h 


H 


H 


J ljtn 

J 0> J ! 

K 


m 

n 

N 

Q 

r 

r 

o 

r 

w 

R 


s 


S 


-Height of instantaneous free surface at tank wall = h J a . 
-Value of h at vapor ingestion. 

-Thickness of viscous film = h Q /a . 

-Mean curvature of free surface = aH . 

-Mean curvature of equilibrium free surface at the container 

axis = aH . 

o 

th - t 

-m zero of . 

-Bessel functions of first kind. 

-Ratio of nondime ns ional draining time to nondimensional 
sloshing time = {w/(l + B)3^^ . 

-Summation index. 

-Exterior normal direction. 

-Normal coordinate of point on free surface meridian in 
local normal- surface coordinates = N/a . 

-Dimensional volumetric drain rate. 

-Radial coordinate = r/a . 

-Drain radius = r /a . 

o' 

-Value of r at intersection of free surface and container 

wall = r /a . 
w 

-Radial coordinate of point on free surface meridian = R/a . 
-Radial coordinate on computer plots = r/a . 

-Reynolds number = pQ/^a(i . 

-Arc length along free surface or a mesh line "parallel" 
to the free surface. 

-Surface coordinate of point on free surface meridian in 
local normal- surface coordinates = S/a . 
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t 

T 

t* 


u 

U 

U' 


v 

c 

V 

w- 

V 

o 

V 

V 

r 

V 

s 


V 

o . 


VOL 

W 

X 

z 

Z 


Greek Alphabet 

9 

© 
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-Time = Qi/na^ ( nondimens ional draining time). 

-Time on computer plots = Qt/na^ . 

-Time = t{(l + B)a/pa^} (nondimensional sloshing time). 
-Time at vapor ingestion. 

-Auxiliary velocity in s-direction for mesh control. 

-Periodic function. 

-Derivative of U with respect to time. 

-Velocity along container axis. 

-Velocity along container wall. 

-Dimensional average bulk speed of liquid = Q/tra 2 • 

-Liquid volume in tank = v/a? . 

-Residual volume measured in hemispherical volumes. 

-Volume of right cylinder of radius 1 below the hemispherical 

bottom of tank with drain of radius r . 

o 

-Initial liquid volume in tank = . 

-Initial tank volume on computer plots = V q /tt . 

2 2 3 

-Weber number = pQ /tt ca . 

-Generic spatial variable = X/a . 

-Axial coordinate = z/a . 

-Axial coordinate of point on free surface meridian = Z/a . 
-Axial coordinate on computer plot = z/a . 


-Azimuthal angle m right cylindrical coordinates. 
-Contact angle (= 5° for numerical calculations). 
-Generic period of wave. 
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-Period of wave = \Q/na^ . 

X* -Period of wave = X( (l + B)cr/pa^)^^ . 

(i -Dimensional liquid viscosity.- 

p -Dimensional liquid density, 

cr -Dimensional surface tension. 

9 -Velocity potential = ua 9 /Q, (nondimensional draining potential). 

-Auxiliary velocity potential. 

9^ -9^ = 9 - 9^ (residual velocity potential). 

- r 

<p* -Velocity potential = 91 (l + B)oa/pj (nondimensional 

sloshing potential). 

9 n -Exterior normal derivative of 9 . 

«0 s -Tangential derivative of 9, positive to left of 9 n . 


Miscellaneous 

V 


-Gradient operator = aV . 
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